Model selection for Ch 1 - eelgrass and sand flat data

*NOTES:*

  1. Decided we should include site as a random effect for all purse seine models because the # of sites is well above 5 and so we deemed it important to acknowledge the nested nature of our sampling design.
  1. Early on in model selection we determined that our data has a negative binomial distribution. All models started with Poisson distributions and all were overdispersed, which was corrected by the neg binomial distribution. All models here use negative binomial distribution with a log link.
  1. We found that the relationship between abundance and Julian day is commonly better described as quadratic (vs. linear), and have incorporated this into the global model (i.e including Julian day as two parameters: Jday + Jday^2). In some cases this term was removed from the global model if it did not converge (see purse.chinook.Rmd)
  1. Habitat variables were tested for collinearity, and through this process we narrowed them down to Leaf Area Index and Mean Turbidity only to include in the global model. These two variables still do covary, but make biological sense as separate variables.

4.5) Site E5 (South Ferry) was missing measurements for leaf area index, so the values from site E6 (Point Roberts) were used. Site E7 (near coal port in intercauseway) was missing values for turbidity, so the mean of sites E3 and E4 (Intercauseway N and S) were used. Archipelago values for leaf area index were used for sites E3, E4, and E7. All eelgrass measurements were set to 0 for sand flat sites. Mean turbidity values were calculated for all sites (including sand flat)

  1. We used AICc to rank all models in dredge, as it penalizes the strength of likelihood more for very small sample sizes, and is more adaptive to moderate sample sizes, so that we could apply the same information criterion across all models (Harrison et al 2018; Brewer et al 2016) – could explain better…
  1. We used model averaging following a delta AICc of less than 4, [[and followed the nesting rule to remove any artificially added models that were represented by smaller models with better AICc scores (Harrison et al 2018).]]- didnt do this because weights were low enough anyway so arbitrary. Kept at delta 4 because weights dropped off considerably before even that point. The choice of cuttoff is somewhat arbitrary, but within the accepted standard of ranges to incorporate enough uncertainty while also having some selectivity. From Jenn:

Hi Lia,

I don’t have a reference that explicitly says that you should use a delta AIC cut off of 4, but my decision to use 4 as a cutoff was mainly based on two things:

  1. In Burnham and Anderson’s 2002 book they say “… an alternative rule of thumb for an approximate 95% confidence set on the K-L best model is the subset of all models having delta less than or equal to some value that is roughly 4-7” (p. 170), which seems to suggest that any value within that range would be a reasonable choice. I also found that I had a lot more models in the top model set when using a higher delta AIC cut off value (e.g. 7), so 4 seemed like a better choice.
  1. A lot of the stats for my paper are modeled after the stats used in another coral reef-related paper by Emily Darling (https://link.springer.com/article/10.1007/s00338-017-1539-z), where they also used a delta AIC cut off of 4 (although they don’t provide a reference for it).

Chinook purse seine model

response = average abundance of Chinook salmon [per purse seine site]
This is 1st of 4 separate purse models for (1)Chinook, (2)chin, (3)other migratory species and (4)resident species

Load all data - standardize variables, create subgroups from purse seine data

# load Purse seine data aggregated by species/set; contains unidentified
# larval
purse <- read.csv("/Users/Lia/Documents/Git/Fraser-salmon/all.data/purse.catch.csv")

# get summary stats for each model parameter for appendix
summary(purse)
##       Year               Date        Temp.surf        DO.surf     
##  Min.   :0.0000   6/2/2017 : 139   Min.   : 7.49   Min.   : 43.1  
##  1st Qu.:0.0000   5/19/2017: 111   1st Qu.:12.51   1st Qu.:101.0  
##  Median :0.0000   7/15/2017:  99   Median :14.06   Median :112.2  
##  Mean   :0.4481   6/3/2017 :  90   Mean   :13.90   Mean   :110.1  
##  3rd Qu.:1.0000   6/20/2017:  87   3rd Qu.:15.31   3rd Qu.:119.1  
##  Max.   :1.0000   24-Jun-16:  81   Max.   :20.41   Max.   :156.5  
##                   (Other)  :1620                                  
##    DOmg.surf        pH.surf         Sal.surf        Temp.bot    
##  Min.   : 4.04   Min.   :7.780   Min.   : 0.55   Min.   : 6.19  
##  1st Qu.: 9.25   1st Qu.:8.060   1st Qu.:12.96   1st Qu.:11.84  
##  Median :10.17   Median :8.250   Median :23.02   Median :13.22  
##  Mean   :10.08   Mean   :8.268   Mean   :20.63   Mean   :13.33  
##  3rd Qu.:10.99   3rd Qu.:8.520   3rd Qu.:28.80   3rd Qu.:14.79  
##  Max.   :14.56   Max.   :9.000   Max.   :33.92   Max.   :20.22  
##                                                                 
##      DO.bot         DOmg.bot         pH.bot        Sal.bot     
##  Min.   : 45.3   Min.   : 4.13   Min.   :7.02   Min.   : 1.31  
##  1st Qu.:100.6   1st Qu.: 9.16   1st Qu.:8.03   1st Qu.:21.11  
##  Median :110.0   Median : 9.79   Median :8.22   Median :27.17  
##  Mean   :111.1   Mean   :10.09   Mean   :8.23   Mean   :24.61  
##  3rd Qu.:121.5   3rd Qu.:11.19   3rd Qu.:8.41   3rd Qu.:29.97  
##  Max.   :177.3   Max.   :16.85   Max.   :9.20   Max.   :34.05  
##                                                                
##      J.date         SIN_TIME           COS_TIME           Round       
##  Min.   : 79.0   Min.   :-0.99771   Min.   :-0.9999   Min.   : 1.000  
##  1st Qu.:133.0   1st Qu.: 0.01075   1st Qu.:-0.9765   1st Qu.: 5.000  
##  Median :154.0   Median : 0.47276   Median :-0.8644   Median : 6.000  
##  Mean   :161.3   Mean   : 0.32371   Mean   :-0.6879   Mean   : 6.173  
##  3rd Qu.:182.0   3rd Qu.: 0.75370   3rd Qu.:-0.5185   3rd Qu.: 8.000  
##  Max.   :286.0   Max.   : 0.99999   Max.   : 0.2102   Max.   :10.000  
##                                                                       
##       Habitat          Site          Set       
##  Eelgrass :1435   E5     :275   Min.   :1.000  
##  Sand flat: 792   E3     :265   1st Qu.:1.000  
##                   E2     :255   Median :2.000  
##                   E4     :235   Mean   :2.007  
##                   E1     :205   3rd Qu.:3.000  
##                   SF6    :178   Max.   :4.000  
##                   (Other):814   NA's   :1      
##                      Species       Life_stage    abundance      
##  Shiner surfperch        :329   0       : 82   Min.   :   0.00  
##  Three-spined stickleback:320   adult   :884   1st Qu.:   1.00  
##  Starry flounder         :293   juvenile:962   Median :   2.00  
##  Bay pipefish            :147   NA's    :299   Mean   :  19.77  
##  Chinook                 :131                  3rd Qu.:   6.00  
##  Surf smelt              : 84                  Max.   :1698.00  
##  (Other)                 :923                                   
##        class          round              month    
##  0        :  82   Min.   : 1.000   April    :357  
##  migratory: 370   1st Qu.: 5.000   July     :341  
##  resident :1717   Median : 6.000   June     :538  
##  NA's     :  58   Mean   : 6.078   March    :153  
##                   3rd Qu.: 7.500   May      :656  
##                   Max.   :10.000   October  : 80  
##                                    September:102
sapply(purse, sd)
##       Year       Date  Temp.surf    DO.surf  DOmg.surf    pH.surf 
##  0.4974146 14.0883494  2.5111848 15.2915703  1.4735729  0.2599412 
##   Sal.surf   Temp.bot     DO.bot   DOmg.bot     pH.bot    Sal.bot 
##  9.4074621  2.4150753 17.5070536  1.6823201  0.2550196  7.5839695 
##     J.date   SIN_TIME   COS_TIME      Round    Habitat       Site 
## 44.7545523  0.5426650  0.3573595  2.2216665  0.4788129  3.8004176 
##        Set    Species Life_stage  abundance      class      round 
##         NA 12.4404013         NA 80.1448604         NA  2.0653150 
##      month 
##  1.6950018
# eelgrass summary
eelgrass <- subset(purse, Habitat %in% "Eelgrass")
summary(eelgrass)
##       Year               Date       Temp.surf        DO.surf     
##  Min.   :0.0000   6/2/2017 :139   Min.   : 7.49   Min.   : 57.5  
##  1st Qu.:0.0000   5/19/2017:111   1st Qu.:12.23   1st Qu.:105.0  
##  Median :0.0000   7/15/2017: 99   Median :13.85   Median :113.6  
##  Mean   :0.4544   6/20/2017: 87   Mean   :13.64   Mean   :111.9  
##  3rd Qu.:1.0000   24-Jun-16: 81   3rd Qu.:15.04   3rd Qu.:122.7  
##  Max.   :1.0000   10-Jul-16: 77   Max.   :18.85   Max.   :156.5  
##                   (Other)  :841                                  
##    DOmg.surf         pH.surf         Sal.surf        Temp.bot    
##  Min.   : 5.530   Min.   :7.780   Min.   :10.38   Min.   : 7.07  
##  1st Qu.: 9.190   1st Qu.:8.110   1st Qu.:22.81   1st Qu.:11.63  
##  Median :10.160   Median :8.300   Median :27.36   Median :13.06  
##  Mean   : 9.965   Mean   :8.313   Mean   :25.56   Mean   :13.06  
##  3rd Qu.:10.880   3rd Qu.:8.550   3rd Qu.:29.64   3rd Qu.:14.65  
##  Max.   :14.190   Max.   :9.000   Max.   :33.92   Max.   :17.53  
##                                                                  
##      DO.bot         DOmg.bot         pH.bot        Sal.bot     
##  Min.   : 64.9   Min.   : 5.96   Min.   :7.02   Min.   :14.89  
##  1st Qu.:102.7   1st Qu.: 9.06   1st Qu.:8.12   1st Qu.:26.76  
##  Median :113.1   Median :10.02   Median :8.24   Median :29.14  
##  Mean   :113.3   Mean   :10.16   Mean   :8.26   Mean   :28.50  
##  3rd Qu.:124.3   3rd Qu.:11.20   3rd Qu.:8.42   3rd Qu.:30.83  
##  Max.   :177.3   Max.   :16.85   Max.   :9.20   Max.   :33.88  
##                                                                
##      J.date         SIN_TIME          COS_TIME           Round       
##  Min.   : 79.0   Min.   :-0.9882   Min.   :-0.9953   Min.   : 1.000  
##  1st Qu.:134.0   1st Qu.:-0.1606   1st Qu.:-0.9765   1st Qu.: 5.000  
##  Median :153.0   Median : 0.4878   Median :-0.8278   Median : 7.000  
##  Mean   :162.7   Mean   : 0.3044   Mean   :-0.6873   Mean   : 6.271  
##  3rd Qu.:192.0   3rd Qu.: 0.7423   3rd Qu.:-0.5037   3rd Qu.: 8.000  
##  Max.   :286.0   Max.   : 1.0000   Max.   : 0.2102   Max.   :10.000  
##                                                                      
##       Habitat          Site          Set       
##  Eelgrass :1435   E5     :275   Min.   :1.000  
##  Sand flat:   0   E3     :265   1st Qu.:1.000  
##                   E2     :255   Median :2.000  
##                   E4     :235   Mean   :2.019  
##                   E1     :205   3rd Qu.:3.000  
##                   E7     :105   Max.   :4.000  
##                   (Other): 95   NA's   :1      
##                      Species       Life_stage    abundance      
##  Three-spined stickleback:245   0       : 17   Min.   :   0.00  
##  Shiner surfperch        :242   adult   :665   1st Qu.:   1.00  
##  Starry flounder         :162   juvenile:539   Median :   2.00  
##  Bay pipefish            :131   NA's    :214   Mean   :  26.05  
##  Tubesnout               : 83                  3rd Qu.:  10.00  
##  Surf smelt              : 74                  Max.   :1698.00  
##  (Other)                 :498                                   
##        class          round              month    
##  0        :  17   Min.   : 1.000   April    :215  
##  migratory: 214   1st Qu.: 5.000   July     :231  
##  resident :1185   Median : 6.000   June     :347  
##  NA's     :  19   Mean   : 6.147   March    : 95  
##                   3rd Qu.: 8.000   May      :418  
##                   Max.   :10.000   October  : 52  
##                                    September: 77
sapply(eelgrass, sd)
##       Year       Date  Temp.surf    DO.surf  DOmg.surf    pH.surf 
##  0.4980858 13.9351400  2.3601026 15.3515577  1.3885806  0.2706839 
##   Sal.surf   Temp.bot     DO.bot   DOmg.bot     pH.bot    Sal.bot 
##  6.0302167  2.2570560 18.3091380  1.7728437  0.2680020  3.4087278 
##     J.date   SIN_TIME   COS_TIME      Round    Habitat       Site 
## 45.5103676  0.5526942  0.3602235  2.2162626  0.0000000  1.7656842 
##        Set    Species Life_stage  abundance      class      round 
##         NA 11.6310130         NA 95.3126364         NA  2.0735395 
##      month 
##  1.7078374
# Sand flat summary
sandflat <- subset(purse, Habitat %in% "Sand flat")
summary(sandflat)
##       Year               Date       Temp.surf        DO.surf     
##  Min.   :0.0000   6/3/2017 : 90   Min.   : 7.55   Min.   : 43.1  
##  1st Qu.:0.0000   7/14/2017: 57   1st Qu.:12.62   1st Qu.: 98.4  
##  Median :0.0000   12-Jun-16: 51   Median :14.21   Median :109.9  
##  Mean   :0.4369   6/19/2017: 50   Mean   :14.38   Mean   :106.8  
##  3rd Qu.:1.0000   5/18/2017: 48   3rd Qu.:15.75   3rd Qu.:116.0  
##  Max.   :1.0000   5/6/2017 : 42   Max.   :20.41   Max.   :143.5  
##                   (Other)  :454                                  
##    DOmg.surf        pH.surf         Sal.surf        Temp.bot    
##  Min.   : 4.04   Min.   :7.810   Min.   : 0.55   Min.   : 6.19  
##  1st Qu.: 9.41   1st Qu.:8.015   1st Qu.: 5.27   1st Qu.:12.35  
##  Median :10.36   Median :8.140   Median :10.97   Median :13.65  
##  Mean   :10.28   Mean   :8.187   Mean   :11.70   Mean   :13.83  
##  3rd Qu.:11.44   3rd Qu.:8.310   3rd Qu.:16.42   3rd Qu.:15.50  
##  Max.   :14.56   Max.   :8.700   Max.   :30.84   Max.   :20.22  
##                                                                 
##      DO.bot         DOmg.bot          pH.bot         Sal.bot     
##  Min.   : 45.3   Min.   : 4.130   Min.   :7.770   Min.   : 1.31  
##  1st Qu.: 97.9   1st Qu.: 9.258   1st Qu.:8.020   1st Qu.:12.17  
##  Median :108.2   Median : 9.740   Median :8.120   Median :16.14  
##  Mean   :107.1   Mean   : 9.968   Mean   :8.174   Mean   :17.55  
##  3rd Qu.:114.7   3rd Qu.:11.140   3rd Qu.:8.280   3rd Qu.:24.35  
##  Max.   :144.2   Max.   :13.360   Max.   :8.670   Max.   :34.05  
##                                                                  
##      J.date         SIN_TIME           COS_TIME           Round       
##  Min.   : 79.0   Min.   :-0.99771   Min.   :-0.9999   Min.   : 1.000  
##  1st Qu.:126.0   1st Qu.: 0.01075   1st Qu.:-0.9696   1st Qu.: 4.000  
##  Median :154.0   Median : 0.47276   Median :-0.8644   Median : 6.000  
##  Mean   :158.7   Mean   : 0.35872   Mean   :-0.6889   Mean   : 5.996  
##  3rd Qu.:182.0   3rd Qu.: 0.82719   3rd Qu.:-0.5185   3rd Qu.: 8.000  
##  Max.   :285.0   Max.   : 0.99999   Max.   : 0.2102   Max.   :10.000  
##                                                                       
##       Habitat         Site          Set       
##  Eelgrass :  0   SF6    :178   Min.   :1.000  
##  Sand flat:792   SF1    :137   1st Qu.:1.000  
##                  SF3    :135   Median :2.000  
##                  SF4    :121   Mean   :1.986  
##                  SF2    :117   3rd Qu.:3.000  
##                  SF5    :104   Max.   :3.000  
##                  (Other):  0                  
##                      Species       Life_stage    abundance      
##  Starry flounder         :131   0       : 65   Min.   :  0.000  
##  Shiner surfperch        : 87   adult   :219   1st Qu.:  1.000  
##  Three-spined stickleback: 75   juvenile:423   Median :  1.000  
##  Chinook                 : 74   NA's    : 85   Mean   :  8.374  
##  0                       : 65                  3rd Qu.:  3.250  
##  Arrow goby              : 64                  Max.   :626.000  
##  (Other)                 :296                                   
##        class         round              month    
##  0        : 65   Min.   : 1.000   April    :142  
##  migratory:156   1st Qu.: 4.000   July     :110  
##  resident :532   Median : 6.000   June     :191  
##  NA's     : 39   Mean   : 5.951   March    : 58  
##                  3rd Qu.: 7.500   May      :238  
##                  Max.   :10.000   October  : 28  
##                                   September: 25
sapply(sandflat, sd)
##       Year       Date  Temp.surf    DO.surf  DOmg.surf    pH.surf 
##  0.4963119 14.3669995  2.7012171 14.6284763  1.5976603  0.2170129 
##   Sal.surf   Temp.bot     DO.bot   DOmg.bot     pH.bot    Sal.bot 
##  7.6917446  2.6065195 15.1503214  1.4979394  0.2192653  7.9716065 
##     J.date   SIN_TIME   COS_TIME      Round    Habitat       Site 
## 43.2564944  0.5225303  0.3523348  2.2218854  0.0000000  1.7858508 
##        Set    Species Life_stage  abundance      class      round 
##  0.8083377 13.6352066         NA 37.4836395         NA  2.0455862 
##      month 
##  1.6705429
### Grab just chinook catch - automatically removes any 0s or unidentified
### species.  1
p.1 <- purse[which(purse$Species == "Chinook"), ]


### standardize variables install.packages('robustHD')
### install.packages('robustbase', repos='http://R-Forge.R-project.org')
library(robustHD)

p.1$s.temp <- standardize(p.1$Temp.surf, centerFun = mean, scaleFun = sd)
summary(p.1$s.temp)  #can see that mean is now 0, and SD is on the original scale of x (temp) -- i.e. predictor centering 
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -2.9418 -0.5635 -0.1433  0.0000  0.4847  2.8797
p.1$s.sal <- standardize(p.1$Sal.surf, centerFun = mean, scaleFun = sd)
p.1$s.do <- standardize(p.1$DOmg.surf, centerFun = mean, scaleFun = sd)
p.1$s.pH <- standardize(p.1$pH.surf, centerFun = mean, scaleFun = sd)
p.1$s.J.date <- standardize(p.1$J.date, centerFun = mean, scaleFun = sd)
### Create variable j2 which is Julian day squared in order to represent J
### date as quadratic relationship instead of linear
p.1$j2 <- p.1$s.J.date^2

# summarize by site-day
library(plyr)
p.chin <- ddply(p.1, .(Year, J.date, s.J.date, j2, Habitat, Site, s.temp, s.sal, 
    s.do, s.pH), summarize, abundance = sum(abundance))

# plot abundance~J.date to see data distribution
plot(abundance ~ J.date, data = p.chin)

Habitat variables

Add eelgrass and sand flat habitat variables for each site

p.hab <- read.csv("/Users/Lia/Documents/Git/Fraser-salmon/all.data/site.char.eelgrass.csv")
summary(p.hab)
##     Collection.date         Site.ID       Site     shoot_den     
##             :7                  :6   E1     :1   Min.   :   0.0  
##  18-Jul-16  :2      017 FR      :1   E2     :1   1st Qu.:   0.0  
##  20-Jul-16  :1      017BFR      :1   E3     :1   Median : 304.0  
##  archipelago:3      018 FR      :1   E4     :1   Mean   : 349.6  
##                     030 FR      :1   E5     :1   3rd Qu.: 508.4  
##                     WP126(19 FR):1   E6     :1   Max.   :1118.9  
##                     (Other)     :2   (Other):7                   
##   percent_cov      leaf_area_index  epiphyte_wt        meanturb     
##  Min.   :  0.000   Min.   :0.000   Min.   : 0.000   Min.   : 2.017  
##  1st Qu.:  0.000   1st Qu.:0.000   1st Qu.: 0.000   1st Qu.: 3.465  
##  Median :  7.667   Median :2.000   Median : 0.000   Median : 8.971  
##  Mean   : 39.863   Mean   :2.017   Mean   : 1.788   Mean   : 8.113  
##  3rd Qu.: 75.000   3rd Qu.:3.451   3rd Qu.: 0.484   3rd Qu.:10.941  
##  Max.   :100.000   Max.   :5.831   Max.   :10.112   Max.   :17.810  
##                                    NA's   :3
sapply(p.hab, sd)
## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm): Calling var(x) on a factor x is deprecated and will become an error.
##   Use something like 'all(duplicated(x)[-1L])' to test for a constant vector.

## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm): Calling var(x) on a factor x is deprecated and will become an error.
##   Use something like 'all(duplicated(x)[-1L])' to test for a constant vector.

## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm): Calling var(x) on a factor x is deprecated and will become an error.
##   Use something like 'all(duplicated(x)[-1L])' to test for a constant vector.
## Collection.date         Site.ID            Site       shoot_den 
##        1.290994        2.577019        3.894440      415.152611 
##     percent_cov leaf_area_index     epiphyte_wt        meanturb 
##       44.244763        2.217159              NA        5.011330
# eelgrass summary
eelgrass <- subset(p.hab, Site %in% c("E1", "E2", "E3", "E4", "E5", "E6", "E7"))
summary(eelgrass)
##     Collection.date         Site.ID       Site     shoot_den     
##             :1      017 FR      :1   E1     :1   Min.   : 304.0  
##  18-Jul-16  :2      017BFR      :1   E2     :1   1st Qu.: 400.0  
##  20-Jul-16  :1      018 FR      :1   E3     :1   Median : 508.4  
##  archipelago:3      030 FR      :1   E4     :1   Mean   : 649.2  
##                     WP126(19 FR):1   E5     :1   3rd Qu.: 906.6  
##                     WP227(20 FR):1   E6     :1   Max.   :1118.9  
##                     (Other)     :1   (Other):1                   
##   percent_cov      leaf_area_index  epiphyte_wt        meanturb    
##  Min.   :  7.667   Min.   :2.000   Min.   : 0.484   Min.   :2.017  
##  1st Qu.: 75.000   1st Qu.:2.815   1st Qu.: 0.484   1st Qu.:2.295  
##  Median : 75.000   Median :3.451   Median : 3.640   Median :3.465  
##  Mean   : 74.032   Mean   :3.746   Mean   : 4.469   Mean   :4.430  
##  3rd Qu.: 92.778   3rd Qu.:4.656   3rd Qu.: 7.625   3rd Qu.:5.820  
##  Max.   :100.000   Max.   :5.831   Max.   :10.112   Max.   :9.295  
##                                    NA's   :3
sapply(eelgrass, sd)
## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm): Calling var(x) on a factor x is deprecated and will become an error.
##   Use something like 'all(duplicated(x)[-1L])' to test for a constant vector.

## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm): Calling var(x) on a factor x is deprecated and will become an error.
##   Use something like 'all(duplicated(x)[-1L])' to test for a constant vector.

## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm): Calling var(x) on a factor x is deprecated and will become an error.
##   Use something like 'all(duplicated(x)[-1L])' to test for a constant vector.
## Collection.date         Site.ID            Site       shoot_den 
##        1.214986        2.160247        2.160247      343.146423 
##     percent_cov leaf_area_index     epiphyte_wt        meanturb 
##       31.049192        1.508139              NA        2.771784
# Sand flat summary
sandflat <- subset(p.hab, Site %in% c("SF1", "SF2", "SF3", "SF4", "SF5", "SF6"))
summary(sandflat)
##     Collection.date         Site.ID       Site     shoot_den  percent_cov
##             :6                  :6   SF1    :1   Min.   :0   Min.   :0   
##  18-Jul-16  :0      017 FR      :0   SF2    :1   1st Qu.:0   1st Qu.:0   
##  20-Jul-16  :0      017BFR      :0   SF3    :1   Median :0   Median :0   
##  archipelago:0      018 FR      :0   SF4    :1   Mean   :0   Mean   :0   
##                     030 FR      :0   SF5    :1   3rd Qu.:0   3rd Qu.:0   
##                     WP126(19 FR):0   SF6    :1   Max.   :0   Max.   :0   
##                     (Other)     :0   (Other):0                           
##  leaf_area_index  epiphyte_wt    meanturb     
##  Min.   :0       Min.   :0    Min.   : 8.971  
##  1st Qu.:0       1st Qu.:0    1st Qu.:10.427  
##  Median :0       Median :0    Median :11.901  
##  Mean   :0       Mean   :0    Mean   :12.411  
##  3rd Qu.:0       3rd Qu.:0    3rd Qu.:13.435  
##  Max.   :0       Max.   :0    Max.   :17.810  
## 
sapply(sandflat, sd)
## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm): Calling var(x) on a factor x is deprecated and will become an error.
##   Use something like 'all(duplicated(x)[-1L])' to test for a constant vector.

## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm): Calling var(x) on a factor x is deprecated and will become an error.
##   Use something like 'all(duplicated(x)[-1L])' to test for a constant vector.

## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm): Calling var(x) on a factor x is deprecated and will become an error.
##   Use something like 'all(duplicated(x)[-1L])' to test for a constant vector.
## Collection.date         Site.ID            Site       shoot_den 
##        0.000000        0.000000        1.870829        0.000000 
##     percent_cov leaf_area_index     epiphyte_wt        meanturb 
##        0.000000        0.000000        0.000000        3.145279
## Now add to p.chin
p.hab2 <- p.hab[, c(3:6, 8)]
p.chin.hab <- p.chin
p.chin.hab$leaf_area_index <- p.hab2[match(p.chin.hab$Site, p.hab2$Site), 4]
p.chin.hab$meanturb <- p.hab2[match(p.chin.hab$Site, p.hab2$Site), 5]

## standardize
p.chin.hab$leaf_area_index <- standardize(p.chin.hab$leaf_area_index, centerFun = mean, 
    scaleFun = sd)
p.chin.hab$leaf_area_index <- as.numeric(p.chin.hab$leaf_area_index)
p.chin.hab$meanturb <- standardize(p.chin.hab$meanturb, centerFun = mean, scaleFun = sd)
p.chin.hab$meanturb <- as.numeric(p.chin.hab$meanturb)

VIF for collinearity of habitat variables

Assess variance inflation factors

library(car)
library(GGally)

vif(lm(abundance ~ Habitat + leaf_area_index + s.J.date + j2 + meanturb + Year + 
    s.temp + s.sal + s.do + s.pH, data = p.chin.hab))
##         Habitat leaf_area_index        s.J.date              j2 
##        7.160505        6.972103        3.050891        2.062968 
##        meanturb            Year          s.temp           s.sal 
##        4.898655        1.358753        3.174750        5.458973 
##            s.do            s.pH 
##        2.720738        1.968218

# alias function identifies covariates that are multiples of each other - in
# this case we are OK.
alias(lm(abundance ~ Habitat + leaf_area_index + s.J.date + j2 + meanturb + 
    Year + s.temp + s.sal + s.do + s.pH, data = p.chin.hab))
## Model :
## abundance ~ Habitat + leaf_area_index + s.J.date + j2 + meanturb + 
##     Year + s.temp + s.sal + s.do + s.pH

## Pearson Corr with all vars
year <- p.chin.hab$Year
jday <- p.chin.hab$s.J.date
j2 <- p.chin.hab$j2
temp <- p.chin.hab$s.temp
do <- p.chin.hab$s.do
ph <- p.chin.hab$s.pH
sal <- p.chin.hab$s.sal
lai <- p.chin.hab$leaf_area_index
turb <- p.chin.hab$meanturb
hab <- p.chin.hab$Habitat

habcovar <- cbind.data.frame(hab, year, lai, turb, jday, j2, temp, do, ph, sal)
ggpairs(data = na.omit(habcovar), title = "Pearson Correlation plot habitat variables")

ggsave("/Users/Lia/Documents/Git/Fraser-salmon/Ch1.Seasonal_diversity/expl.plots/PearsonCorrAllVars_chinook.pdf")

Model selection

Full model

We started with a full global model with all the abiotic factors and all the habitat covariates included as well as site as a random effect. That model did not converge, so we performed manual backward selection beginning with the variables with highest VIF, and retaining the top variables from the corresponding marsh model, until the model converged. Throughout the process we viewed the VIF and model assumption plots, and used a hypothesis-based approach to select variables (as well as trying multiple iterations).

library(MASS)
library(lme4)
## Warning: package 'lme4' was built under R version 3.4.4
## Chinook: full model with habitat variables
p.chin1 <- glmer.nb(abundance ~ Habitat + s.J.date + s.temp + Year + s.do + 
    s.pH + (1 | Site), data = p.chin.hab, na.action = "na.fail")
summary(p.chin1)  ## AIC 365.7 
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(4.4747)  ( log )
## Formula: abundance ~ Habitat + s.J.date + s.temp + Year + s.do + s.pH +  
##     (1 | Site)
##    Data: p.chin.hab
## 
##      AIC      BIC   logLik deviance df.resid 
##    365.7    387.8   -173.9    347.7       77 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.4021 -0.5944 -0.2815  0.4967  3.1490 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.1119   0.3345  
## Number of obs: 86, groups:  Site, 13
## 
## Fixed effects:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       0.98589    0.23431   4.208 2.58e-05 ***
## HabitatSand flat -0.51129    0.29608  -1.727  0.08419 .  
## s.J.date         -0.12703    0.13666  -0.930  0.35261    
## s.temp            0.07661    0.12192   0.628  0.52977    
## Year              0.52378    0.18808   2.785  0.00535 ** 
## s.do              0.26861    0.11149   2.409  0.01598 *  
## s.pH             -0.03519    0.12091  -0.291  0.77100    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) HbttSf s.J.dt s.temp Year   s.do  
## HabttSndflt -0.680                                   
## s.J.date    -0.034  0.068                            
## s.temp       0.185 -0.286 -0.604                     
## Year        -0.533  0.084 -0.057  0.121              
## s.do         0.027 -0.112  0.459 -0.084  0.103       
## s.pH        -0.299  0.467 -0.009 -0.331 -0.051 -0.182
## Introduce function to calculate Variance Inflation Factors for GLMMs
vif.lme <- function(fit) {
    ## adapted from rms::vif
    v <- vcov(fit)
    nam <- names(fixef(fit))
    ## exclude intercepts
    ns <- sum(1 * (nam == "Intercept" | nam == "(Intercept)"))
    if (ns > 0) {
        v <- v[-(1:ns), -(1:ns), drop = FALSE]
        nam <- nam[-(1:ns)]
    }
    d <- diag(v)^0.5
    v <- diag(solve(v/(d %o% d)))
    names(v) <- nam
    v
}

# test VIF of full model -- mostly OK
vif.lme(p.chin1)
## HabitatSand flat         s.J.date           s.temp             Year 
##         1.345749         2.229278         2.053339         1.053817 
##             s.do             s.pH 
##         1.409247         1.449275
# Test fit and assumptions of full (global) model Plot residuals vs. fitted
# values
plot(fitted(p.chin1), resid(p.chin1), main = "Global Chinook GLMM", xlab = "Fitted", 
    ylab = "Pearson residuals")

## q plot
qqnorm(x = p.chin.hab$abundance, y = resid(p.chin1), main = "Q-Q Global Chinook GLMM")
qqline(resid(p.chin1), col = 2)

library(sjPlot)
## Warning: package 'sjPlot' was built under R version 3.4.3
## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.2.12
## Current Matrix version is 1.2.11
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
# QQ plot of random effects quantiles against standard normal quantiles
sjp.glmer(p.chin1, type = "re.qq")

# Check model assumptions - note that we are not actually assuming linear
# relationships with Jday, but can't tell this general function that.
sjp.glmer(p.chin1, type = "ma")

# See marginal and conditional pseudo R2 for full model
library(piecewiseSEM)
rsquared(p.chin1, aicc = TRUE)
##      Class            Family Link  n  Marginal Conditional
## 1 glmerMod Negative Binomial  log 86 0.3587852   0.5289181

Model selection - dredge

Then continued to dredge function to determine optimal model using AICc (chin has somewhat small sample size: 86 observations at 13 sites).

library(MuMIn)
library(knitr)
# Generate model set
p.model.set.chin <- dredge(p.chin1)
# p.model.set.chin.AIC<- dredge(p.chin1, rank = 'AIC') #compare to AIC to
# see if much different

# Create model selection table
p.model_table.chin <- model.sel(p.model.set.chin)
options(scipen = 7)
names(p.model_table.chin) <- c("(Int)", "Habitat", "DO", "Jday", "pH", "Temp", 
    "Year", "df", "LL", "AICc", "delta", "weight")
kable(head(p.model_table.chin, n = 100), digits = 3)
(Int) Habitat DO Jday pH Temp Year df LL AICc delta weight
36 0.956 + 0.313 NA NA NA 0.514 6 -174.344 361.751 0.000 0.20456373180
35 0.705 NA 0.308 NA NA NA 0.544 5 -175.782 362.314 0.563 0.15440510347
40 0.953 + 0.273 -0.078 NA NA 0.510 7 -174.060 363.556 1.805 0.08297937494
39 0.705 NA 0.262 -0.087 NA NA 0.539 6 -175.433 363.929 2.177 0.06886675942
44 0.972 + 0.316 NA -0.032 NA 0.514 7 -174.302 364.039 2.288 0.06516961517
52 0.956 + 0.312 NA NA -0.003 0.513 7 -174.343 364.123 2.371 0.06250457753
51 0.705 NA 0.302 NA NA -0.026 0.539 6 -175.739 364.542 2.790 0.05068528339
43 0.705 NA 0.307 NA 0.019 NA 0.544 6 -175.768 364.599 2.847 0.04926704491
56 0.965 + 0.263 -0.127 NA 0.065 0.521 8 -173.899 365.669 3.917 0.02885127572
48 0.958 + 0.275 -0.076 -0.010 NA 0.510 8 -174.056 365.982 4.231 0.02466810618
47 0.705 NA 0.253 -0.098 0.045 NA 0.539 7 -175.354 366.145 4.393 0.02274471869
55 0.705 NA 0.257 -0.113 NA 0.034 0.545 7 -175.388 366.212 4.461 0.02199037633
60 0.977 + 0.318 NA -0.036 0.009 0.516 8 -174.297 366.465 4.713 0.01938115146
59 0.705 NA 0.297 NA 0.037 -0.037 0.536 7 -175.691 366.818 5.067 0.01624125323
4 1.310 + 0.288 NA NA NA NA 5 -178.186 367.122 5.371 0.01395158938
37 0.726 NA NA -0.213 NA NA 0.489 5 -178.409 367.569 5.817 0.01115858447
38 0.955 + NA -0.208 NA NA 0.461 6 -177.270 367.604 5.852 0.01096561749
64 0.986 + 0.269 -0.127 -0.035 0.077 0.524 9 -173.857 368.082 6.331 0.00863202856
3 1.022 NA 0.281 NA NA NA NA 4 -179.981 368.457 6.705 0.00715789797
63 0.705 NA 0.251 -0.113 0.038 0.022 0.543 8 -175.337 368.545 6.793 0.00685048720
8 1.305 + 0.242 -0.088 NA NA NA 6 -177.819 368.701 6.950 0.00633389601
54 0.976 + NA -0.281 NA 0.109 0.482 7 -176.789 369.013 7.262 0.00541882828
45 0.725 NA NA -0.226 0.098 NA 0.490 6 -178.010 369.084 7.333 0.00523102988
20 1.301 + 0.283 NA NA -0.027 NA 6 -178.140 369.344 7.592 0.00459433372
53 0.725 NA NA -0.266 NA 0.080 0.505 6 -178.142 369.346 7.595 0.00458792518
12 1.326 + 0.292 NA -0.031 NA NA 6 -178.144 369.351 7.599 0.00457835700
46 0.927 + NA -0.216 0.055 NA 0.461 7 -177.147 369.730 7.978 0.00378748302
7 1.020 NA 0.230 -0.095 NA NA NA 5 -179.557 369.864 8.112 0.00354241144
19 1.017 NA 0.271 NA NA -0.047 NA 5 -179.840 370.429 8.678 0.00266946629
11 1.022 NA 0.279 NA 0.013 NA NA 5 -179.974 370.698 8.947 0.00233361241
24 1.314 + 0.237 -0.111 NA 0.032 NA 7 -177.779 370.994 9.243 0.00201258153
16 1.308 + 0.244 -0.086 -0.007 NA NA 7 -177.817 371.070 9.318 0.00193809015
34 0.965 + NA NA NA NA 0.457 5 -180.173 371.096 9.344 0.00191330723
61 0.724 NA NA -0.258 0.080 0.052 0.501 7 -177.910 371.256 9.505 0.00176582713
33 0.731 NA NA NA NA NA 0.482 4 -181.446 371.385 9.634 0.00165528238
62 0.966 + NA -0.280 0.017 0.102 0.481 8 -176.779 371.427 9.676 0.00162075377
6 1.275 + NA -0.205 NA NA NA 5 -180.404 371.557 9.806 0.00151883333
28 1.314 + 0.287 NA -0.021 -0.019 NA 7 -178.124 371.684 9.933 0.00142544201
15 1.020 NA 0.220 -0.105 0.042 NA NA 6 -179.488 372.039 10.288 0.00119358174
23 1.021 NA 0.229 -0.098 NA 0.004 NA 6 -179.556 372.176 10.424 0.00111490008
5 1.013 NA NA -0.208 NA NA NA 4 -181.879 372.251 10.499 0.00107376227
27 1.015 NA 0.264 NA 0.046 -0.063 NA 6 -179.765 372.594 10.842 0.00090454585
49 0.728 NA NA NA NA -0.078 0.471 5 -181.003 372.756 11.004 0.00083425616
50 0.949 + NA NA NA -0.061 0.448 6 -179.912 372.887 11.136 0.00078114591
22 1.300 + NA -0.253 NA 0.076 NA 6 -180.158 373.379 11.628 0.00061073715
41 0.732 NA NA NA 0.054 NA 0.482 5 -181.318 373.386 11.635 0.00060864901
42 0.960 + NA NA 0.011 NA 0.457 6 -180.167 373.398 11.647 0.00060498902
32 1.327 + 0.240 -0.111 -0.021 0.039 NA 8 -177.764 373.399 11.648 0.00060478186
14 1.250 + NA -0.212 0.052 NA NA 6 -180.287 373.637 11.886 0.00053689177
13 1.014 NA NA -0.220 0.093 NA NA 5 -181.519 373.789 12.037 0.00049766863
57 0.730 NA NA NA 0.105 -0.109 0.464 6 -180.569 374.200 12.449 0.00040510387
21 1.019 NA NA -0.238 NA 0.049 NA 5 -181.774 374.298 12.547 0.00038577707
31 1.019 NA 0.222 -0.098 0.046 -0.012 NA 7 -179.483 374.401 12.650 0.00036636411
2 1.278 + NA NA NA NA NA 4 -183.286 375.066 13.314 0.00026280648
58 0.921 + NA NA 0.050 -0.077 0.445 7 -179.818 375.072 13.320 0.00026202510
30 1.284 + NA -0.251 0.027 0.066 NA 7 -180.131 375.699 13.947 0.00019152803
1 1.013 NA NA NA NA NA NA 3 -184.837 375.966 14.214 0.00016757094
29 1.016 NA NA -0.229 0.086 0.016 NA 6 -181.510 376.082 14.331 0.00015809118
18 1.253 + NA NA NA -0.075 NA 5 -182.893 376.536 14.784 0.00012602500
17 1.003 NA NA NA NA -0.092 NA 4 -184.246 376.986 15.235 0.00010059887
10 1.273 + NA NA 0.011 NA NA 5 -183.281 377.312 15.560 0.00008549940
9 1.013 NA NA NA 0.049 NA NA 4 -184.730 377.953 16.201 0.00006204732
25 1.000 NA NA NA 0.113 -0.127 NA 5 -183.753 378.256 16.504 0.00005332741
26 1.218 + NA NA 0.060 -0.094 NA 6 -182.760 378.583 16.831 0.00004528767

# determine at which delta AICc score you reach cumulative AICc weight of
# 0.95 (Harrison et al. 2018)
find_cumsum = function(df, delta) {
    for (i in nrow(df)) {
        df = df[df$delta >= delta, ]
        return(min(df$delta[cumsum(df$weight) >= 0.95]))
    }
}
find_cumsum(p.model_table.chin, 0)  ##THIS SUGGESTS AVERAGING UP TO delta AICc 7.59. We will use a cutoff of delta AICc < 4.
## [1] 7.592112

# Model averaging Version 1: use all models with delta AIC score less than 4
p.model.set.chin.4 <- get.models(p.model.set.chin, subset = delta < 4)
p.avg_model.chin.4 <- model.avg(p.model.set.chin.4)
summary(p.avg_model.chin.4)
## 
## Call:
## model.avg(object = p.model.set.chin.4)
## 
## Component model call: 
## glmer(formula = abundance ~ <9 unique rhs>, data = p.chin.hab, 
##      family = negative.binomial(theta = 4.47467009362264), na.action = 
##      na.fail)
## 
## Component models: 
##       df  logLik   AICc delta weight
## 126    6 -174.34 361.75  0.00   0.27
## 26     5 -175.78 362.31  0.56   0.20
## 1236   7 -174.06 363.56  1.80   0.11
## 236    6 -175.43 363.93  2.18   0.09
## 1246   7 -174.30 364.04  2.29   0.08
## 1256   7 -174.34 364.12  2.37   0.08
## 256    6 -175.74 364.54  2.79   0.07
## 246    6 -175.77 364.60  2.85   0.06
## 12356  8 -173.90 365.67  3.92   0.04
## 
## Term codes: 
##  Habitat     s.do s.J.date     s.pH   s.temp     Year 
##        1        2        3        4        5        6 
## 
## Model-averaged coefficients:  
## (full average) 
##                    Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)       0.8516625  0.2406532   0.2433611   3.500 0.000466 ***
## HabitatSand flat -0.2651547  0.2999038   0.3019090   0.878 0.379802    
## s.do              0.3002695  0.0993223   0.1007970   2.979 0.002892 ** 
## Year              0.5255259  0.1874715   0.1903397   2.761 0.005763 ** 
## s.J.date         -0.0210789  0.0659719   0.0666429   0.316 0.751778    
## s.pH             -0.0015055  0.0434506   0.0440854   0.034 0.972758    
## s.temp            0.0004921  0.0428127   0.0434058   0.011 0.990954    
##  
## (conditional average) 
##                   Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)       0.851662   0.240653    0.243361   3.500 0.000466 ***
## HabitatSand flat -0.458153   0.258816    0.262814   1.743 0.081288 .  
## s.do              0.300269   0.099322    0.100797   2.979 0.002892 ** 
## Year              0.525526   0.187472    0.190340   2.761 0.005763 ** 
## s.J.date         -0.089507   0.111159    0.112846   0.793 0.427672    
## s.pH             -0.010094   0.112125    0.113774   0.089 0.929303    
## s.temp            0.002658   0.099476    0.100855   0.026 0.978971    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Relative variable importance: 
##                      s.do Year Habitat s.J.date s.temp s.pH
## Importance:          1.00 1.00 0.58    0.24     0.19   0.15
## N containing models:    9    9    5       3        3      2
p.chin.ci <- data.frame(confint(p.avg_model.chin.4, full = TRUE))

# Get pseudo R squared values for models up to delta < 4 model.list.chin<-
# list(#manually list top x models from dredge)
p.chin.4.Rsq <- rsquared(p.model.set.chin.4, aicc = TRUE)

## write tables to .csv for easy comparison and plugging into appendix table
p.avg_model_4df.chin <- data.frame(p.avg_model.chin.4$msTable)
p.avg_model_components4.chin <- cbind(p.chin.4.Rsq, p.avg_model_4df.chin)
p.r = data.frame(Coeff = rownames(p.avg_model_4df.chin, rep(NA, length(p.avg_model_components4.chin))))
p.avg_model_components4.chin <- cbind(p.avg_model_components4.chin, p.r)
p.avg_model_components4.chin <- p.avg_model_components4.chin[, -c(7, 8)]

# write.csv(p.avg_model_components4.chin,
# '/Users/Lia/Documents/Git/Fraser-salmon/all.data/avg_model_components4_pursechin.csv')

Parameter Plot

The results of model averaging including all top ranked models up to delta AICc 4

## Warning: package 'cowplot' was built under R version 3.4.3
## 
## Attaching package: 'cowplot'
## The following objects are masked from 'package:sjPlot':
## 
##     plot_grid, save_plot
## The following object is masked from 'package:ggplot2':
## 
##     ggsave

Check fit of top model by avg – including all parameters with weight >0.5

# Test fit and assumptions of final (averaged) model, following the nesting
# rule

p.avg.chin <- glmer.nb(abundance ~ s.do + Year + Habitat + (1 | Site), data = p.chin.hab, 
    na.action = "na.fail")

summary(p.avg.chin)  # AIC 360.7
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(4.3106)  ( log )
## Formula: abundance ~ s.do + Year + Habitat + (1 | Site)
##    Data: p.chin.hab
## 
##      AIC      BIC   logLik deviance df.resid 
##    360.7    375.4   -174.3    348.7       80 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3717 -0.6216 -0.2597  0.3799  3.1145 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.1021   0.3195  
## Number of obs: 86, groups:  Site, 13
## 
## Fixed effects:
##                  Estimate Std. Error z value  Pr(>|z|)    
## (Intercept)       0.95601    0.21976   4.350 0.0000136 ***
## s.do              0.31274    0.09416   3.322  0.000895 ***
## Year              0.51413    0.18717   2.747  0.006017 ** 
## HabitatSand flat -0.45215    0.25286  -1.788  0.073757 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) s.do   Year  
## s.do        -0.032              
## Year        -0.602  0.112       
## HabttSndflt -0.634 -0.047  0.148
vif.lme(p.avg.chin)
##             s.do             Year HabitatSand flat 
##         1.016995         1.037407         1.026673
### Plot residuals vs. fitted values
plot(fitted(p.avg.chin), resid(p.avg.chin), main = "Averaged Chinook GLMM", 
    xlab = "Fitted", ylab = "Pearson residuals")

## q plot
qqnorm(x = p.chin.hab$abundance, y = resid(p.avg.chin), main = "Q-Q Averaged Chinook GLMM")
qqline(resid(p.avg.chin), col = 2)

library(sjPlot)
# QQ plot of random effects quantiles against standard normal quantiles
sjp.glmer(p.avg.chin, type = "re.qq")

Chum purse seine model

response = average abundance of chum salmon [per purse seine site]
This is 2nd of 4 purse models for (1)Chinook, (2)chum, (3)other migratory species and (4)resident species

Load all data - standardize variables, create subgroups from purse seine data

### Grab just chum catch - automatically removes any 0s or unidentified
### species.  2
p.2 <- purse[which(purse$Species == "Chum"), ]


### standardize variables
p.2$s.temp <- standardize(p.2$Temp.surf, centerFun = mean, scaleFun = sd)
summary(p.2$s.temp)  #can see that mean is now 0, and SD is on the original scale of x (temp) -- i.e. predictor centering 
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -1.9135 -0.1349  0.1520  0.0000  0.4476  2.3575
p.2$s.sal <- standardize(p.2$Sal.surf, centerFun = mean, scaleFun = sd)
p.2$s.do <- standardize(p.2$DOmg.surf, centerFun = mean, scaleFun = sd)
p.2$s.pH <- standardize(p.2$pH.surf, centerFun = mean, scaleFun = sd)
p.2$s.J.date <- standardize(p.2$J.date, centerFun = mean, scaleFun = sd)
### Create variable j2 which is Julian day squared in order to represent J
### date as quadratic relationship instead of linear
p.2$j2 <- p.2$s.J.date^2

# summarize by site-day
p.chum <- ddply(p.2, .(Year, J.date, s.J.date, j2, Habitat, Site, s.temp, s.sal, 
    s.do, s.pH), summarize, abundance = sum(abundance))

# plot abundance~J.date to see data distribution
plot(abundance ~ J.date, data = p.chum)

Habitat variables

Add eelgrass and sand flat habitat variables for each site

## Add habitat variables to p.chum
p.chum.hab <- p.chum
p.chum.hab$leaf_area_index <- p.hab2[match(p.chum.hab$Site, p.hab2$Site), 4]
p.chum.hab$meanturb <- p.hab2[match(p.chum.hab$Site, p.hab2$Site), 5]

## standardize
p.chum.hab$leaf_area_index <- standardize(p.chum.hab$leaf_area_index, centerFun = mean, 
    scaleFun = sd)
p.chum.hab$leaf_area_index <- as.numeric(p.chum.hab$leaf_area_index)
p.chum.hab$meanturb <- standardize(p.chum.hab$meanturb, centerFun = mean, scaleFun = sd)
p.chum.hab$meanturb <- as.numeric(p.chum.hab$meanturb)

VIF for collinearity of habitat variables

Assess variance inflation factors

vif(lm(abundance ~ Habitat + leaf_area_index + s.J.date + j2 + meanturb + Year + 
    s.temp + s.sal + s.do + s.pH, data = p.chum.hab))
##         Habitat leaf_area_index        s.J.date              j2 
##        5.869989        7.362660        3.996986        2.553199 
##        meanturb            Year          s.temp           s.sal 
##        6.961885        1.318973        5.240322        3.951466 
##            s.do            s.pH 
##        2.674625        3.008126

# alias function identifies covariates that are multiples of each other - in
# this case we are OK.
alias(lm(abundance ~ Habitat + leaf_area_index + s.J.date + j2 + meanturb + 
    Year + s.temp + s.sal + s.do + s.pH, data = p.chum.hab))
## Model :
## abundance ~ Habitat + leaf_area_index + s.J.date + j2 + meanturb + 
##     Year + s.temp + s.sal + s.do + s.pH

## Pearson Corr with all vars
year <- p.chum.hab$Year
jday <- p.chum.hab$s.J.date
j2 <- p.chum.hab$j2
temp <- p.chum.hab$s.temp
do <- p.chum.hab$s.do
ph <- p.chum.hab$s.pH
sal <- p.chum.hab$s.sal
lai <- p.chum.hab$leaf_area_index
turb <- p.chum.hab$meanturb
hab <- p.chum.hab$Habitat

habcovar <- cbind.data.frame(hab, year, lai, turb, jday, j2, temp, do, ph, sal)
ggpairs(data = na.omit(habcovar), title = "Pearson Correlation plot habitat variables")

ggsave("/Users/Lia/Documents/Git/Fraser-salmon/Ch1.Seasonal_diversity/expl.plots/PearsonCorrAllVars_chum.pdf")

Model selection

Full model

We started with a full global model with all the abiotic factors and all the habitat covariates included as well as site as a random effect. That model did not converge, so we performed manual backward selection beginning with the variables with highest VIF, and retaining the top variables from the corresponding marsh model, until the model converged. Throughout the process we viewed the VIF and model assumption plots, and used a hypothesis-based approach to select variables (as well as trying multiple iterations).

## Chum: full model with habitat variables
p.chum.1 <- glmer.nb(abundance ~ Habitat + s.J.date + j2 + Year + s.do + meanturb + 
    (1 | Site), data = p.chum.hab, na.action = "na.fail")
summary(p.chum.1)  ## AIC 217.9
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(1.733)  ( log )
## Formula: abundance ~ Habitat + s.J.date + j2 + Year + s.do + meanturb +  
##     (1 | Site)
##    Data: p.chum.hab
## 
##      AIC      BIC   logLik deviance df.resid 
##    217.9    230.5   -100.0    199.9       21 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.1241 -0.6075 -0.1823  0.4802  2.3779 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 1.547    1.244   
## Number of obs: 30, groups:  Site, 12
## 
## Fixed effects:
##                  Estimate Std. Error z value  Pr(>|z|)    
## (Intercept)       1.08269    0.94034   1.151   0.24957    
## HabitatSand flat  0.04549    1.39273   0.033   0.97394    
## s.J.date         -0.98209    0.31976  -3.071   0.00213 ** 
## j2               -0.91609    0.23062  -3.972 0.0000712 ***
## Year              1.75525    0.63867   2.748   0.00599 ** 
## s.do             -0.32625    0.26171  -1.247   0.21253    
## meanturb         -0.87781    0.65204  -1.346   0.17822    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) HbttSf s.J.dt j2     Year   s.do  
## HabttSndflt -0.668                                   
## s.J.date    -0.035  0.026                            
## j2          -0.128 -0.016  0.528                     
## Year        -0.546 -0.033 -0.097 -0.088              
## s.do        -0.073 -0.186 -0.106  0.463  0.222       
## meanturb     0.528 -0.780  0.207  0.134 -0.041  0.090
# test VIF of full model -- several parameters are a bit high
vif.lme(p.chum.1)
## HabitatSand flat         s.J.date               j2             Year 
##         2.915721         1.914494         2.374107         1.120226 
##             s.do         meanturb 
##         1.874700         2.967015
# Test fit and assumptions of full (global) model Plot residuals vs. fitted
# values
plot(fitted(p.chum.1), resid(p.chum.1), main = "Global Chum1 GLMM", xlab = "Fitted", 
    ylab = "Pearson residuals")

## q plot
qqnorm(x = p.chum.hab$abundance, y = resid(p.chum.1), main = "Q-Q Global Chum GLMM")
qqline(resid(p.chum.1), col = 2)

# QQ plot of random effects quantiles against standard normal quantiles --
# this is odd, all at 0???
sjp.glmer(p.chum.1, type = "re.qq")

# Check model assumptions
sjp.glmer(p.chum.1, type = "ma")

# See marginal and conditional pseudo R2 for full model - quite high.
rsquared(p.chum.1, aicc = TRUE)
##      Class            Family Link  n  Marginal Conditional
## 1 glmerMod Negative Binomial  log 30 0.4332422   0.9569078
## Note did an iteration of global model that also converged: habitat + Jday
## + Year + DO + meanturb + Temp, AIC was higher (225.1), VIFs were higher
## and temp relationship was non linear, but marginal r squared was over 0.9

Model selection - dredge

Then continued to dredge function to determine optimal model using AICc (chum has smaller sample size: 30 observations at 12 sites).

# Generate model set p.model.set.chum <- dredge(p.chum.1, extra =
# 'r.squaredGLMM') #NOTE r.squaredGLMM doesn't work with lme4. Could
# re-model with lme...
p.model.set.chum <- dredge(p.chum.1)

# Create model selection table
p.model_table.chum <- model.sel(p.model.set.chum)
options(scipen = 7)
names(p.model_table.chum) <- c("(Int)", "Habitat", "Jday^2", "Mean turbidity", 
    "Dissolved oxygen", "Jday", "Year", "df", "LL", "AICc", "delta", "weight")
kable(head(p.model_table.chum, n = 100), digits = 3)
(Int) Habitat Jday^2 Mean turbidity Dissolved oxygen Jday Year df LL AICc delta weight
55 0.872 NA -0.780 -0.908 NA -1.031 1.929 7 -100.763 220.617 0.000 0.27914960100
51 0.876 NA -0.695 NA NA -0.906 1.741 6 -103.281 222.215 1.598 0.12553014162
52 1.658 + -0.767 NA NA -0.975 1.852 7 -101.577 222.244 1.627 0.12372240894
63 1.103 NA -0.916 -0.861 -0.324 -0.983 1.757 8 -99.963 222.784 2.167 0.09444822990
59 1.165 NA -0.881 NA -0.418 -0.841 1.554 7 -101.971 223.032 2.415 0.08343453825
56 1.009 + -0.783 -0.802 NA -1.029 1.929 8 -100.738 224.333 3.717 0.04352769299
60 1.785 + -0.899 NA -0.324 -0.930 1.682 8 -100.838 224.533 3.917 0.03938599388
27 2.522 NA -0.857 NA -0.533 -0.798 NA 6 -105.076 225.804 5.187 0.02087050651
43 1.130 NA -0.528 NA -0.458 NA 1.445 6 -105.359 226.370 5.754 0.01571996457
23 2.494 NA -0.675 -0.813 NA -0.986 NA 6 -105.381 226.414 5.797 0.01538258219
35 0.803 NA -0.297 NA NA NA 1.672 5 -107.047 226.594 5.978 0.01405462651
31 2.614 NA -0.883 -0.753 -0.470 -0.910 NA 7 -103.769 226.629 6.012 0.01381356356
19 2.367 NA -0.608 NA NA -0.880 NA 5 -107.144 226.788 6.172 0.01275418758
64 1.082 + -0.916 -0.878 -0.326 -0.982 1.756 9 -99.963 226.926 6.309 0.01190731187
33 0.714 NA NA NA NA NA 1.435 4 -108.752 227.104 6.487 0.01089345745
39 0.784 NA -0.319 -0.490 NA NA 1.800 6 -105.854 227.361 6.744 0.00958002895
20 3.140 + -0.667 NA NA -0.944 NA 6 -105.933 227.519 6.902 0.00885228460
28 3.124 + -0.871 NA -0.472 -0.869 NA 7 -104.398 227.887 7.271 0.00736241268
36 1.256 + -0.335 NA NA NA 1.763 6 -106.129 227.910 7.293 0.00728029986
37 0.700 NA NA -0.454 NA NA 1.515 5 -107.848 228.195 7.578 0.00631279883
47 1.077 NA -0.508 -0.433 -0.399 NA 1.549 7 -104.627 228.344 7.728 0.00585855514
11 2.379 NA -0.518 NA -0.575 NA NA 5 -108.085 228.670 8.053 0.00497902565
34 1.056 + NA NA NA NA 1.470 5 -108.254 229.009 8.392 0.00420310884
44 1.409 + -0.514 NA -0.399 NA 1.514 7 -104.992 229.076 8.459 0.00406387588
49 0.691 NA NA NA NA -0.195 1.402 5 -108.465 229.429 8.813 0.00340557955
24 2.610 + -0.678 -0.723 NA -0.986 NA 7 -105.366 229.823 9.206 0.00279745513
41 0.738 NA NA NA -0.051 NA 1.396 5 -108.719 229.937 9.321 0.00264174320
1 1.997 NA NA NA NA NA NA 3 -111.646 230.216 9.599 0.00229830952
32 2.510 + -0.883 -0.836 -0.476 -0.908 NA 8 -103.758 230.374 9.757 0.00212353075
53 0.673 NA NA -0.529 NA -0.243 1.479 6 -107.398 230.448 9.831 0.00204684565
40 0.922 + -0.327 -0.387 NA NA 1.803 7 -105.819 230.730 10.113 0.00177759393
3 2.210 NA -0.217 NA NA NA NA 4 -110.666 230.931 10.315 0.00160712724
15 2.405 NA -0.506 -0.372 -0.545 NA NA 6 -107.656 230.964 10.347 0.00158144855
38 0.655 + NA -0.490 NA NA 1.516 6 -107.844 231.340 10.723 0.00131012190
45 0.702 NA NA -0.454 -0.004 NA 1.512 6 -107.847 231.347 10.730 0.00130567347
12 2.625 + -0.510 NA -0.547 NA NA 6 -107.921 231.493 10.877 0.00121340915
50 1.066 + NA NA NA -0.208 1.435 6 -107.926 231.505 10.888 0.00120636769
5 2.026 NA NA -0.405 NA NA NA 4 -111.037 231.673 11.057 0.00110890044
17 1.940 NA NA NA NA -0.228 NA 4 -111.225 232.050 11.433 0.00091868531
48 0.995 + -0.509 -0.500 -0.407 NA 1.544 8 -104.616 232.089 11.473 0.00090071322
42 1.055 + NA NA 0.007 NA 1.477 6 -108.254 232.160 11.543 0.00086962621
2 2.324 + NA NA NA NA NA 4 -111.299 232.199 11.582 0.00085265679
9 1.958 NA NA NA -0.153 NA NA 4 -111.346 232.291 11.675 0.00081411494
7 2.265 NA -0.228 -0.425 NA NA NA 5 -109.938 232.376 11.759 0.00078039039
57 0.675 NA NA NA 0.033 -0.215 1.422 6 -108.454 232.561 11.944 0.00071147652
4 2.638 + -0.240 NA NA NA NA 5 -110.115 232.731 12.114 0.00065356950
21 1.976 NA NA -0.486 NA -0.269 NA 5 -110.448 233.395 12.779 0.00046881044
61 0.604 NA NA -0.571 0.129 -0.325 1.569 7 -107.253 233.597 12.981 0.00042374699
54 0.572 + NA -0.610 NA -0.247 1.481 7 -107.382 233.855 13.239 0.00037247648
13 2.002 NA NA -0.390 -0.126 NA NA 5 -110.824 234.148 13.531 0.00032183655
18 2.304 + NA NA NA -0.245 NA 5 -110.827 234.155 13.538 0.00032066042
16 2.227 + -0.507 -0.514 -0.556 NA NA 7 -107.619 234.329 13.713 0.00029386438
6 2.000 + NA -0.431 NA NA NA 5 -111.035 234.569 13.953 0.00026063609
58 1.062 + NA NA 0.129 -0.287 1.523 7 -107.792 234.675 14.058 0.00024726358
10 2.246 + NA NA -0.119 NA NA 5 -111.118 234.737 14.120 0.00023974300
46 0.652 + NA -0.494 -0.009 NA 1.508 7 -107.843 234.777 14.160 0.00023498507
25 1.935 NA NA NA -0.076 -0.176 NA 5 -111.170 234.840 14.223 0.00022769712
8 2.365 + -0.233 -0.350 NA NA NA 6 -109.924 235.501 14.884 0.00016362014
22 1.881 + NA -0.563 NA -0.272 NA 6 -110.436 236.524 15.907 0.00009810346
29 1.975 NA NA -0.481 -0.014 -0.259 NA 6 -110.446 236.544 15.927 0.00009711378
14 1.861 + NA -0.500 -0.137 NA NA 6 -110.799 237.250 16.633 0.00006822990
26 2.297 + NA NA -0.017 -0.229 NA 6 -110.825 237.302 16.685 0.00006648867
62 0.579 + NA -0.592 0.126 -0.325 1.568 8 -107.252 237.361 16.745 0.00006452894
30 1.865 + NA -0.567 -0.023 -0.256 NA 7 -110.431 239.953 19.337 0.00001765871

# Model averaging Version 1: use all models with delta AIC score less than 4
p.model.set.chum.4 <- get.models(p.model.set.chum, subset = delta < 4)
p.avg_model.chum.4 <- model.avg(p.model.set.chum.4)
summary(p.avg_model.chum.4)
## 
## Call:
## model.avg(object = p.model.set.chum.4)
## 
## Component model call: 
## glmer(formula = abundance ~ <7 unique rhs>, data = p.chum.hab, 
##      family = negative.binomial(theta = 1.73301192961737), na.action = 
##      na.fail)
## 
## Component models: 
##       df  logLik   AICc delta weight
## 2356   7 -100.76 220.62  0.00   0.35
## 256    6 -103.28 222.22  1.60   0.16
## 1256   7 -101.58 222.24  1.63   0.16
## 23456  8  -99.96 222.78  2.17   0.12
## 2456   7 -101.97 223.03  2.42   0.11
## 12356  8 -100.74 224.33  3.72   0.06
## 12456  8 -100.84 224.53  3.92   0.05
## 
## Term codes: 
##  Habitat       j2 meanturb     s.do s.J.date     Year 
##        1        2        3        4        5        6 
## 
## Model-averaged coefficients:  
## (full average) 
##                  Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)       1.10770    0.79318     0.83059   1.334 0.182324    
## j2               -0.79740    0.22508     0.23645   3.372 0.000745 ***
## meanturb         -0.46845    0.53746     0.54738   0.856 0.392108    
## s.J.date         -0.97137    0.32964     0.34732   2.797 0.005162 ** 
## Year              1.81439    0.64083     0.67526   2.687 0.007210 ** 
## HabitatSand flat -0.34304    0.80815     0.82582   0.415 0.677851    
## s.do             -0.09921    0.21202     0.21707   0.457 0.647652    
##  
## (conditional average) 
##                  Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)        1.1077     0.7932      0.8306   1.334 0.182324    
## j2                -0.7974     0.2251      0.2364   3.372 0.000745 ***
## meanturb          -0.8863     0.4197      0.4433   1.999 0.045590 *  
## s.J.date          -0.9714     0.3296      0.3473   2.797 0.005162 ** 
## Year               1.8144     0.6408      0.6753   2.687 0.007210 ** 
## HabitatSand flat  -1.3102     1.1078      1.1565   1.133 0.257274    
## s.do              -0.3604     0.2630      0.2776   1.298 0.194214    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Relative variable importance: 
##                      j2   s.J.date Year meanturb s.do Habitat
## Importance:          1.00 1.00     1.00 0.53     0.28 0.26   
## N containing models:    7    7        7    3        3    3
p.chum.ci <- data.frame(confint(p.avg_model.chum.4, full = TRUE))

# Get pseudo R squared values for models up to delta < 4
p.chum.4.Rsq <- rsquared(p.model.set.chum.4, aicc = TRUE)

## write tables to .csv for easy comparison and plugging into appendix table
p.avg_model_4df.chum <- data.frame(p.avg_model.chum.4$msTable)
p.avg_model_components4.chum <- cbind(p.chum.4.Rsq, p.avg_model_4df.chum)
p.r = data.frame(Coeff = rownames(p.avg_model_4df.chum, rep(NA, length(p.avg_model_components4.chum))))
p.avg_model_components4.chum <- cbind(p.avg_model_components4.chum, p.r)
p.avg_model_components4.chum <- p.avg_model_components4.chum[, -c(7, 8)]
# write.csv(p.avg_model_components4.chum,
# '/Users/Lia/Documents/Git/Fraser-salmon/all.data/avg_model_components4_pursechum.csv')

Parameter Plot

The results of model averaging including all top ranked models up to delta AICc 4

Migratory purse seine model

response = average abundance of migra salmon [per purse seine site]
This is 3rd of 4 purse models for (1)Chinook, (2)migra, (3)other migratory species and (4)resident species

Load all data - standardize variables, create subgroups from purse seine data

### Grab just migratory catch - automatically removes any 0s or unidentified
### species.  3
p.3 <- purse[which(purse$class == "migratory"), ]
p.3 <- p.3[which(p.3$Species != "Chinook"), ]
p.3 <- p.3[which(p.3$Species != "Chum"), ]

### standardize variables

p.3$s.temp <- standardize(p.3$Temp.surf, centerFun = mean, scaleFun = sd)
summary(p.3$s.temp)  #can see that mean is now 0, and SD is on the original scale of x (temp) -- i.e. predictor centering 
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -3.49165 -0.69003 -0.05758  0.00000  0.38148  3.14128
p.3$s.sal <- standardize(p.3$Sal.surf, centerFun = mean, scaleFun = sd)
p.3$s.do <- standardize(p.3$DOmg.surf, centerFun = mean, scaleFun = sd)
p.3$s.pH <- standardize(p.3$pH.surf, centerFun = mean, scaleFun = sd)
p.3$s.J.date <- standardize(p.3$J.date, centerFun = mean, scaleFun = sd)
### Create variable j2 which is Julian day squared in order to represent J
### date as quadratic relationship instead of linear
p.3$j2 <- p.3$s.J.date^2

# summarize by site-day
p.migra <- ddply(p.3, .(Year, J.date, s.J.date, j2, Habitat, Site, s.temp, s.sal, 
    s.do, s.pH), summarize, abundance = sum(abundance))

# plot abundance~J.date to see data distribution
plot(abundance ~ J.date, data = p.migra)

Habitat variables

Add eelgrass and sand flat habitat variables for each site

## Add habitat variables to p.migra
p.migra.hab <- p.migra
p.migra.hab$leaf_area_index <- p.hab2[match(p.migra.hab$Site, p.hab2$Site), 
    4]
p.migra.hab$meanturb <- p.hab2[match(p.migra.hab$Site, p.hab2$Site), 5]

## standardize
p.migra.hab$leaf_area_index <- standardize(p.migra.hab$leaf_area_index, centerFun = mean, 
    scaleFun = sd)
p.migra.hab$leaf_area_index <- as.numeric(p.migra.hab$leaf_area_index)
p.migra.hab$meanturb <- standardize(p.migra.hab$meanturb, centerFun = mean, 
    scaleFun = sd)
p.migra.hab$meanturb <- as.numeric(p.migra.hab$meanturb)

VIF for collinearity of habitat variables

Assess variance inflation factors

vif(lm(abundance ~ Habitat + leaf_area_index + s.J.date + j2 + meanturb + Year + 
    s.temp + s.sal + s.do + s.pH, data = p.migra.hab))
##         Habitat leaf_area_index        s.J.date              j2 
##        5.538994        5.839467        4.459692        4.894707 
##        meanturb            Year          s.temp           s.sal 
##        4.645115        1.710438        3.176537        4.043365 
##            s.do            s.pH 
##        2.208030        2.635699

# alias function identifies covariates that are multiples of each other - in
# this case we are OK.
alias(lm(abundance ~ Habitat + leaf_area_index + s.J.date + j2 + meanturb + 
    Year + s.temp + s.sal + s.do + s.pH, data = p.migra.hab))
## Model :
## abundance ~ Habitat + leaf_area_index + s.J.date + j2 + meanturb + 
##     Year + s.temp + s.sal + s.do + s.pH

## Pearson Corr with all vars
year <- p.migra.hab$Year
jday <- p.migra.hab$s.J.date
j2 <- p.migra.hab$j2
temp <- p.migra.hab$s.temp
do <- p.migra.hab$s.do
ph <- p.migra.hab$s.pH
sal <- p.migra.hab$s.sal
lai <- p.migra.hab$leaf_area_index
turb <- p.migra.hab$meanturb
hab <- p.migra.hab$Habitat

habcovar <- cbind.data.frame(hab, year, lai, turb, jday, j2, temp, do, ph, sal)
ggpairs(data = na.omit(habcovar), title = "Pearson Correlation plot habitat variables")

ggsave("/Users/Lia/Documents/Git/Fraser-salmon/Ch1.Seasonal_diversity/expl.plots/PearsonCorrAllVars_migratory.pdf")

Model selection

Full model

We started with a full global model with all the abiotic factors and all the habitat covariates included as well as site as a random effect. That model did not converge, so we performed manual backward selection beginning with the variables with highest VIF, and retaining the top variables from the corresponding marsh model, until the model converged. Throughout the process we viewed the VIF and model assumption plots, and used a hypothesis-based approach to select variables (as well as trying multiple iterations). In this case, the marsh model did include j^2, but looking at the purse migratory abundance distribution over time (plot above), this didn’t make sense. We removed it first and re-ran VIF and found that they all improved, supporting that this was not a good fit for this data. The marsh model also contained temp, but the distribution was quite odd with this subset and the model converged when it was removed so we left it out.

## Migratory: full model with habitat variables
p.migra.1 <- glmer.nb(abundance ~ Habitat + s.J.date + Year + s.do + s.pH + 
    (1 | Site), data = p.migra.hab, na.action = "na.fail")
summary(p.migra.1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(0.6208)  ( log )
## Formula: abundance ~ Habitat + s.J.date + Year + s.do + s.pH + (1 | Site)
##    Data: p.migra.hab
## 
##      AIC      BIC   logLik deviance df.resid 
##    722.8    742.8   -353.4    706.8       82 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.7721 -0.6739 -0.4504  0.1326  4.9788 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.08809  0.2968  
## Number of obs: 90, groups:  Site, 13
## 
## Fixed effects:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        2.6473     0.2626  10.080  < 2e-16 ***
## HabitatSand flat   0.5550     0.3853   1.440  0.14976    
## s.J.date           0.2969     0.1557   1.907  0.05654 .  
## Year               0.2358     0.3525   0.669  0.50360    
## s.do               0.5340     0.1701   3.139  0.00169 ** 
## s.pH              -0.0393     0.1737  -0.226  0.82105    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) HbttSf s.J.dt Year   s.do  
## HabttSndflt -0.440                            
## s.J.date    -0.157 -0.226                     
## Year        -0.562 -0.114  0.348              
## s.do         0.116 -0.322  0.452  0.070       
## s.pH         0.047  0.401 -0.134 -0.335 -0.391
# test VIF of full model
vif.lme(p.migra.1)
## HabitatSand flat         s.J.date             Year             s.do 
##         1.261972         1.521183         1.348178         1.576720 
##             s.pH 
##         1.541817
# Test fit and assumptions of full (global) model Plot residuals vs. fitted
# values
plot(fitted(p.migra.1), resid(p.migra.1), main = "Global Migratory GLMM", xlab = "Fitted", 
    ylab = "Pearson residuals")

## q plot
qqnorm(x = p.migra.hab$abundance, y = resid(p.migra.1), main = "Q-Q Global Migratory GLMM")
qqline(resid(p.migra.1), col = 2)

# QQ plot of random effects quantiles against standard normal quantiles
sjp.glmer(p.migra.1, type = "re.qq")

# Check model assumptions
sjp.glmer(p.migra.1, type = "ma")

# See marginal and conditional pseudo R2 for full model
rsquared(p.migra.1, aicc = TRUE)
##      Class            Family Link  n  Marginal Conditional
## 1 glmerMod Negative Binomial  log 90 0.7175173   0.9011521

Model selection - dredge

Then continued to dredge function to determine optimal model using AICc (migra has 90 observations at 13 sites).

# Generate model set
p.model.set.migra <- dredge(p.migra.1)

# Create model selection table
p.model_table.migra <- model.sel(p.model.set.migra)
options(scipen = 7)
names(p.model_table.migra) <- c("(Int)", "Habitat", "Dissolved oxygen", "Jday", 
    "pH", "Year", "df", "LL", "AICc", "delta", "weight")
kable(head(p.model_table.migra, n = 100), digits = 3)
(Int) Habitat Dissolved oxygen Jday pH Year df LL AICc delta weight
7 2.967 NA 0.573 0.308 NA NA 5 -354.746 720.205 0.000 0.1818628771
8 2.749 + 0.526 0.261 NA NA 6 -353.625 720.262 0.057 0.1767502619
4 2.742 + 0.377 NA NA NA 5 -355.279 721.271 1.066 0.1067219127
15 2.961 NA 0.600 0.304 -0.090 NA 6 -354.571 722.154 1.949 0.0686322157
24 2.651 + 0.519 0.292 NA 0.209 7 -353.432 722.229 2.024 0.0661099024
23 2.889 NA 0.570 0.334 NA 0.165 6 -354.629 722.269 2.064 0.0648033627
3 3.015 NA 0.392 NA NA NA 4 -357.053 722.577 2.372 0.0555458697
16 2.749 + 0.526 0.261 0.001 NA 7 -353.625 722.616 2.411 0.0544808307
12 2.749 + 0.386 NA -0.022 NA 6 -355.268 723.549 3.343 0.0341780859
20 2.758 + 0.382 NA NA -0.032 6 -355.273 723.558 3.353 0.0340130156
31 2.835 NA 0.607 0.342 -0.129 0.261 7 -354.302 723.971 3.765 0.0276772373
11 3.012 NA 0.435 NA -0.122 NA 5 -356.719 724.151 3.946 0.0252861080
32 2.647 + 0.534 0.297 -0.039 0.236 8 -353.406 724.590 4.385 0.0203034131
19 3.067 NA 0.406 NA NA -0.111 5 -356.993 724.701 4.495 0.0192126748
2 2.718 + NA NA NA NA 4 -358.536 725.543 5.337 0.0126112635
28 2.757 + 0.387 NA -0.019 -0.017 7 -355.267 725.900 5.694 0.0105487279
27 3.019 NA 0.436 NA -0.120 -0.014 6 -356.718 726.447 6.242 0.0080220936
10 2.675 + NA NA 0.129 NA 5 -358.185 727.083 6.878 0.0058372244
1 3.025 NA NA NA NA NA 3 -360.610 727.499 7.294 0.0047419292
18 2.645 + NA NA NA 0.148 5 -358.416 727.546 7.340 0.0046320854
6 2.719 + NA 0.046 NA NA 5 -358.476 727.666 7.461 0.0043612742
14 2.667 + NA 0.090 0.166 NA 6 -357.972 728.956 8.751 0.0022883439
26 2.656 + NA NA 0.118 0.046 6 -358.175 729.363 9.158 0.0018672219
22 2.603 + NA 0.093 NA 0.243 6 -358.216 729.444 9.238 0.0017933045
5 3.020 NA NA 0.061 NA NA 4 -360.517 729.504 9.298 0.0017402399
17 2.981 NA NA NA NA 0.092 4 -360.568 729.606 9.400 0.0016538281
9 3.025 NA NA NA 0.015 NA 4 -360.605 729.680 9.475 0.0015931678
30 2.607 + NA 0.111 0.139 0.146 7 -357.892 731.149 10.944 0.0007644781
21 2.935 NA NA 0.092 NA 0.175 5 -360.387 731.488 11.283 0.0006451873
13 3.019 NA NA 0.067 0.032 NA 5 -360.496 731.706 11.501 0.0005785645
25 2.980 NA NA NA -0.002 0.094 5 -360.567 731.849 11.644 0.0005386564
29 2.936 NA NA 0.092 0.006 0.171 6 -360.386 733.785 13.579 0.0002046418

# Model averaging Version 1: use all models with delta AIC score less than 4
p.model.set.migra.4 <- get.models(p.model.set.migra, subset = delta < 4)
p.avg_model.migra.4 <- model.avg(p.model.set.migra.4)
summary(p.avg_model.migra.4)
## 
## Call:
## model.avg(object = p.model.set.migra.4)
## 
## Component model call: 
## glmer(formula = abundance ~ <12 unique rhs>, data = p.migra.hab, 
##      family = negative.binomial(theta = 0.620809109139615), na.action = 
##      na.fail)
## 
## Component models: 
##      df  logLik   AICc delta weight
## 23    5 -354.75 720.21  0.00   0.20
## 123   6 -353.63 720.26  0.06   0.20
## 12    5 -355.28 721.27  1.07   0.12
## 234   6 -354.57 722.15  1.95   0.08
## 1235  7 -353.43 722.23  2.02   0.07
## 235   6 -354.63 722.27  2.06   0.07
## 2     4 -357.05 722.58  2.37   0.06
## 1234  7 -353.63 722.62  2.41   0.06
## 124   6 -355.27 723.55  3.34   0.04
## 125   6 -355.27 723.56  3.35   0.04
## 2345  7 -354.30 723.97  3.77   0.03
## 24    5 -356.72 724.15  3.95   0.03
## 
## Term codes: 
##  Habitat     s.do s.J.date     s.pH     Year 
##        1        2        3        4        5 
## 
## Model-averaged coefficients:  
## (full average) 
##                  Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)       2.83840    0.24781     0.25062  11.326  < 2e-16 ***
## s.do              0.50708    0.17048     0.17243   2.941  0.00327 ** 
## s.J.date          0.20944    0.18325     0.18449   1.135  0.25628    
## HabitatSand flat  0.33699    0.41559     0.41800   0.806  0.42013    
## s.pH             -0.01516    0.08198     0.08295   0.183  0.85502    
## Year              0.03416    0.17355     0.17558   0.195  0.84576    
##  
## (conditional average) 
##                  Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)       2.83840    0.24781     0.25062  11.326  < 2e-16 ***
## s.do              0.50708    0.17048     0.17243   2.941  0.00327 ** 
## s.J.date          0.29309    0.14992     0.15204   1.928  0.05388 .  
## HabitatSand flat  0.63942    0.36652     0.37169   1.720  0.08538 .  
## s.pH             -0.06459    0.15954     0.16165   0.400  0.68945    
## Year              0.15890    0.34684     0.35157   0.452  0.65128    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Relative variable importance: 
##                      s.do s.J.date Habitat s.pH Year
## Importance:          1.00 0.71     0.53    0.23 0.21
## N containing models:   12    7        6       5    4
p.migra.ci <- data.frame(confint(p.avg_model.migra.4, full = TRUE))

# Get pseudo R squared values for models up to delta < 4 model.list.migra<-
# list(#manually list top x models from dredge)
p.migra.4.Rsq <- rsquared(p.model.set.migra.4, aicc = TRUE)

## write tables to .csv for easy comparison and plugging into appendix table
p.avg_model_4df.migra <- data.frame(p.avg_model.migra.4$msTable)
p.avg_model_components4.migra <- cbind(p.migra.4.Rsq, p.avg_model_4df.migra)
p.r = data.frame(Coeff = rownames(p.avg_model_4df.migra, rep(NA, length(p.avg_model_components4.migra))))
p.avg_model_components4.migra <- cbind(p.avg_model_components4.migra, p.r)
p.avg_model_components4.migra <- p.avg_model_components4.migra[, -c(7, 8)]
# write.csv(p.avg_model_components4.migra,
# '/Users/Lia/Documents/Git/Fraser-salmon/all.data/avg_model_components4_pursemigratory.csv')

Parameter Plot

The results of model averaging including all top ranked models up to delta AICc 4

Resident purse seine model

response = average abundance of res salmon [per purse seine site]
This is 4th of 4 purse models for (1)Chinook, (2)res, (3)other migratory species and (4)resident species

Load all data - standardize variables, create subgroups from purse seine data

### Grab just resident catch - automatically removes any 0s or unidentified
### species.  4
p.4 <- purse[which(purse$class == "resident"), ]

### standardize variables
p.4$s.temp <- standardize(p.4$Temp.surf, centerFun = mean, scaleFun = sd)
summary(p.4$s.temp)  #can see that mean is now 0, and SD is on the original scale of x (temp) -- i.e. predictor centering 
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -2.64937 -0.61349  0.01961  0.00000  0.55450  2.59038
p.4$s.sal <- standardize(p.4$Sal.surf, centerFun = mean, scaleFun = sd)
p.4$s.do <- standardize(p.4$DOmg.surf, centerFun = mean, scaleFun = sd)
p.4$s.pH <- standardize(p.4$pH.surf, centerFun = mean, scaleFun = sd)
p.4$s.J.date <- standardize(p.4$J.date, centerFun = mean, scaleFun = sd)
### Create variable j2 which is Julian day squared in order to represent J
### date as quadratic relationship instead of linear
p.4$j2 <- p.4$s.J.date^2

# summarize by site-day
p.res <- ddply(p.4, .(Year, J.date, s.J.date, j2, Habitat, Site, s.temp, s.sal, 
    s.do, s.pH), summarize, abundance = sum(abundance))

# plot abundance~J.date to see data distribution
plot(abundance ~ J.date, data = p.res)

Habitat variables

Add eelgrass and sand flat habitat variables for each site

## Add habitat variables to p.res
p.res.hab <- p.res
p.res.hab$leaf_area_index <- p.hab2[match(p.res.hab$Site, p.hab2$Site), 4]
p.res.hab$meanturb <- p.hab2[match(p.res.hab$Site, p.hab2$Site), 5]

## standardize
p.res.hab$leaf_area_index <- standardize(p.res.hab$leaf_area_index, centerFun = mean, 
    scaleFun = sd)
p.res.hab$leaf_area_index <- as.numeric(p.res.hab$leaf_area_index)
p.res.hab$meanturb <- standardize(p.res.hab$meanturb, centerFun = mean, scaleFun = sd)
p.res.hab$meanturb <- as.numeric(p.res.hab$meanturb)

VIF for collinearity of habitat variables

Assess variance inflation factors

vif(lm(abundance ~ Habitat + leaf_area_index + s.J.date + j2 + meanturb + Year + 
    s.temp + s.sal + s.do + s.pH, data = p.res.hab))
##         Habitat leaf_area_index        s.J.date              j2 
##        6.141154        6.308545        3.432246        4.060365 
##        meanturb            Year          s.temp           s.sal 
##        4.315397        1.452055        5.963063        4.107177 
##            s.do            s.pH 
##        1.935367        2.201143

# alias function identifies covariates that are multiples of each other - in
# this case we are OK.
alias(lm(abundance ~ Habitat + leaf_area_index + s.J.date + j2 + meanturb + 
    Year + s.temp + s.sal + s.do + s.pH, data = p.res.hab))
## Model :
## abundance ~ Habitat + leaf_area_index + s.J.date + j2 + meanturb + 
##     Year + s.temp + s.sal + s.do + s.pH

## Pearson Corr with all vars
year <- p.res.hab$Year
jday <- p.res.hab$s.J.date
j2 <- p.res.hab$j2
temp <- p.res.hab$s.temp
do <- p.res.hab$s.do
ph <- p.res.hab$s.pH
sal <- p.res.hab$s.sal
lai <- p.res.hab$leaf_area_index
turb <- p.res.hab$meanturb
hab <- p.res.hab$Habitat

habcovar <- cbind.data.frame(hab, year, lai, turb, jday, j2, temp, do, ph, sal)
ggpairs(data = na.omit(habcovar), title = "Pearson Correlation plot habitat variables")

ggsave("/Users/Lia/Documents/Git/Fraser-salmon/Ch1.Seasonal_diversity/expl.plots/PearsonCorrAllVars_residents.pdf")

Model selection

Full model

We started with a full global model with all the abiotic factors and all the habitat covariates included as well as site as a random effect. Model did not converge - found J2 was biggest contributor (even just Year + Jday + Habitat + j2 did not converge). Once we removed that and the highly correlated habitat vars, model converged.

## Resident: full model with habitat variables
p.res.1 <- glmer.nb(abundance ~ Year + s.J.date + Habitat + s.sal + s.temp + 
    s.do + (1 | Site), data = p.res.hab, na.action = "na.fail")

summary(p.res.1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(0.662)  ( log )
## Formula: abundance ~ Year + s.J.date + Habitat + s.sal + s.temp + s.do +  
##     (1 | Site)
##    Data: p.res.hab
## 
##      AIC      BIC   logLik deviance df.resid 
##   2119.4   2148.5  -1050.7   2101.4      178 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.8062 -0.6744 -0.4913  0.1393  8.9313 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.1537   0.392   
## Number of obs: 187, groups:  Site, 13
## 
## Fixed effects:
##                  Estimate Std. Error z value      Pr(>|z|)    
## (Intercept)        5.9729     0.2325  25.692       < 2e-16 ***
## Year              -0.1753     0.2244  -0.781         0.435    
## s.J.date           0.2195     0.1254   1.751         0.080 .  
## HabitatSand flat  -2.2840     0.3808  -5.998 0.00000000199 ***
## s.sal             -0.1503     0.1695  -0.887         0.375    
## s.temp             0.7147     0.1395   5.125 0.00000029745 ***
## s.do               0.1235     0.1160   1.065         0.287    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Year   s.J.dt HbttSf s.sal  s.temp
## Year        -0.375                                   
## s.J.date    -0.007  0.152                            
## HabttSndflt -0.664  0.024 -0.046                     
## s.sal       -0.386  0.235 -0.126  0.632              
## s.temp      -0.150  0.332 -0.458  0.193  0.429       
## s.do        -0.047  0.053  0.374  0.122  0.237 -0.220
# test VIF of full model
vif.lme(p.res.1)
##             Year         s.J.date HabitatSand flat            s.sal 
##         1.332541         1.601720         1.736491         2.321479 
##           s.temp             s.do 
##         1.904814         1.361182
# Test fit and assumptions of full (global) model Plot residuals vs. fitted
# values
plot(fitted(p.res.1), resid(p.res.1), main = "Global Resident GLMM", xlab = "Fitted", 
    ylab = "Pearson residuals")

## q plot
qqnorm(x = p.res.hab$abundance, y = resid(p.res.1), main = "Q-Q Global Resident GLMM")
qqline(resid(p.res.1), col = 2)

# QQ plot of random effects quantiles against standard normal quantiles
sjp.glmer(p.res.1, type = "re.qq")

# Check model assumptions
sjp.glmer(p.res.1, type = "ma")

# See marginal and conditional pseudo R2 for full model -- very high
rsquared(p.res.1, aicc = TRUE)
##      Class            Family Link   n  Marginal Conditional
## 1 glmerMod Negative Binomial  log 187 0.9187988   0.9962039

Model selection - dredge

Then continued to dredge function to determine optimal model using AICc (probably unnecessary but have done for all others: res has 187 observations at 13 sites).

# Generate model set
p.model.set.res <- dredge(p.res.1)

# Create model selection table
p.model_table.res <- model.sel(p.model.set.res)
options(scipen = 7)
names(p.model_table.res) <- c("(Int)", "Habitat", "Dissolved oxygen", "Jday", 
    "Salinity", "Temp", "Year", "df", "LL", "AICc", "delta", "weight")
kable(head(p.model_table.res, n = 100), digits = 3)
(Int) Habitat Dissolved oxygen Jday Salinity Temp Year df LL AICc delta weight
18 5.863 + NA NA NA 0.908 NA 5 -1053.126 2116.584 0.000 1.531220e-01
22 5.848 + NA 0.148 NA 0.848 NA 6 -1052.116 2116.699 0.115 1.445893e-01
24 5.857 + 0.148 0.218 NA 0.787 NA 7 -1051.291 2117.207 0.622 1.121679e-01
30 5.918 + NA 0.181 -0.161 0.788 NA 7 -1051.618 2117.861 1.277 8.087970e-02
50 5.917 + NA NA NA 0.870 -0.188 6 -1052.740 2117.946 1.362 7.748999e-02
26 5.901 + NA NA -0.078 0.885 NA 6 -1052.996 2118.458 1.874 5.999424e-02
20 5.869 + 0.045 NA NA 0.899 NA 6 -1053.030 2118.526 1.942 5.798267e-02
54 5.884 + NA 0.136 NA 0.827 -0.123 7 -1051.957 2118.540 1.956 5.758843e-02
32 5.908 + 0.128 0.233 -0.120 0.752 NA 8 -1051.026 2118.860 2.276 4.907511e-02
56 5.895 + 0.148 0.206 NA 0.767 -0.127 8 -1051.119 2119.047 2.463 4.469397e-02
62 5.987 + NA 0.168 -0.191 0.746 -0.185 8 -1051.284 2119.378 2.793 3.788309e-02
58 5.997 + NA NA -0.131 0.820 -0.247 7 -1052.401 2119.428 2.843 3.694746e-02
52 5.929 + 0.060 NA NA 0.855 -0.206 7 -1052.574 2119.773 3.189 3.108460e-02
64 5.973 + 0.124 0.220 -0.150 0.715 -0.175 9 -1050.722 2120.462 3.877 2.203178e-02
28 5.897 + 0.031 NA -0.062 0.884 NA 7 -1052.956 2120.537 3.953 2.121642e-02
60 5.994 + 0.037 NA -0.114 0.818 -0.251 8 -1052.341 2121.492 4.908 1.316368e-02
17 4.848 NA NA NA NA 0.909 NA 4 -1063.257 2134.734 18.149 1.753690e-05
21 4.850 NA NA 0.138 NA 0.855 NA 5 -1062.396 2135.123 18.539 1.443497e-05
49 4.928 NA NA NA NA 0.876 -0.191 5 -1062.888 2136.108 19.524 8.819499e-06
23 4.845 NA 0.103 0.188 NA 0.809 NA 6 -1062.012 2136.490 19.906 7.287296e-06
25 4.857 NA NA NA 0.035 0.919 NA 5 -1063.235 2136.801 20.217 6.237934e-06
19 4.847 NA 0.011 NA NA 0.907 NA 5 -1063.251 2136.833 20.249 6.138807e-06
53 4.905 NA NA 0.126 NA 0.836 -0.132 6 -1062.226 2136.919 20.335 5.880206e-06
29 4.838 NA NA 0.148 -0.046 0.838 NA 6 -1062.361 2137.188 20.604 5.139473e-06
51 4.929 NA 0.023 NA NA 0.869 -0.197 6 -1062.864 2138.194 21.610 3.108799e-06
57 4.927 NA NA NA -0.012 0.871 -0.196 6 -1062.886 2138.238 21.654 3.040038e-06
55 4.900 NA 0.102 0.176 NA 0.791 -0.132 7 -1061.841 2138.309 21.724 2.935337e-06
31 4.843 NA 0.101 0.190 -0.009 0.807 NA 7 -1062.010 2138.646 22.062 2.479189e-06
27 4.859 NA 0.023 NA 0.049 0.918 NA 6 -1063.214 2138.894 22.309 2.190767e-06
61 4.896 NA NA 0.140 -0.077 0.804 -0.157 7 -1062.137 2138.899 22.315 2.184566e-06
59 4.930 NA 0.023 NA 0.001 0.870 -0.197 7 -1062.864 2140.353 23.768 1.056271e-06
63 4.896 NA 0.096 0.180 -0.040 0.778 -0.145 8 -1061.819 2140.447 23.863 1.007740e-06
48 6.232 + 0.216 0.608 -0.555 NA -0.470 8 -1064.164 2145.136 28.552 9.661352e-08
46 6.291 + NA 0.509 -0.665 NA -0.525 7 -1065.749 2146.124 29.540 5.895956e-08
16 6.058 + 0.266 0.691 -0.511 NA NA 7 -1066.732 2148.090 31.506 2.206155e-08
14 6.111 + NA 0.587 -0.654 NA NA 6 -1068.989 2150.445 33.861 6.796668e-09
40 5.961 + 0.382 0.679 NA NA -0.357 7 -1069.673 2153.972 37.388 1.165206e-09
8 5.833 + 0.418 0.726 NA NA NA 6 -1071.259 2154.984 38.400 7.025299e-10
42 6.396 + NA NA -0.632 NA -0.746 6 -1074.538 2161.543 44.959 2.644528e-11
38 5.972 + NA 0.499 NA NA -0.444 6 -1074.977 2162.421 45.837 1.704984e-11
44 6.405 + -0.048 NA -0.656 NA -0.749 7 -1074.450 2163.526 46.941 9.814420e-12
6 5.804 + NA 0.544 NA NA NA 5 -1077.488 2165.307 48.723 4.027507e-12
47 4.965 NA 0.212 0.620 -0.499 NA -0.467 7 -1076.309 2167.243 50.659 1.529650e-12
45 4.962 NA NA 0.525 -0.619 NA -0.526 6 -1077.765 2167.996 51.412 1.049658e-12
15 4.787 NA 0.267 0.699 -0.428 NA NA 6 -1078.618 2169.704 53.119 4.469977e-13
13 4.752 NA NA 0.599 -0.585 NA NA 5 -1080.788 2171.907 55.322 1.485640e-13
39 5.037 NA 0.366 0.680 NA NA -0.331 6 -1080.071 2172.608 56.024 1.046052e-13
10 6.097 + NA NA -0.533 NA NA 5 -1081.143 2172.618 56.033 1.041219e-13
7 4.896 NA 0.395 0.723 NA NA NA 5 -1081.323 2172.978 56.394 8.695593e-14
12 6.103 + -0.029 NA -0.551 NA NA 6 -1081.115 2174.697 58.113 3.680686e-14
34 6.030 + NA NA NA NA -0.567 5 -1082.221 2174.773 58.189 3.544188e-14
36 6.033 + 0.114 NA NA NA -0.561 6 -1081.721 2175.908 59.324 2.009358e-14
37 5.077 NA NA 0.523 NA NA -0.410 5 -1084.888 2180.108 63.523 2.460938e-15
2 5.824 + NA NA NA NA NA 4 -1086.409 2181.037 64.453 1.546059e-15
4 5.834 + 0.129 NA NA NA NA 5 -1085.779 2181.889 65.305 1.009782e-15
5 4.903 NA NA 0.566 NA NA NA 4 -1086.866 2181.952 65.368 9.784846e-16
41 5.058 NA NA NA -0.579 NA -0.745 5 -1086.620 2183.571 66.986 4.356292e-16
43 5.052 NA -0.055 NA -0.611 NA -0.749 6 -1086.511 2185.488 68.904 1.670318e-16
33 5.122 NA NA NA NA NA -0.559 4 -1092.449 2193.118 76.534 3.680100e-18
9 4.771 NA NA NA -0.455 NA NA 4 -1092.671 2193.562 76.978 2.948248e-18
35 5.119 NA 0.111 NA NA NA -0.554 5 -1091.977 2194.285 77.701 2.053822e-18
11 4.768 NA -0.022 NA -0.471 NA NA 5 -1092.655 2195.641 79.057 1.042326e-18
1 4.884 NA NA NA NA NA NA 3 -1096.177 2198.486 81.902 2.513613e-19
3 4.883 NA 0.121 NA NA NA NA 4 -1095.621 2199.462 82.878 1.542993e-19

# Model averaging Version 1: use all models with delta AIC score less than 4
p.model.set.res.4 <- get.models(p.model.set.res, subset = delta < 4)
p.avg_model.res.4 <- model.avg(p.model.set.res.4)
summary(p.avg_model.res.4)
## 
## Call:
## model.avg(object = p.model.set.res.4)
## 
## Component model call: 
## glmer(formula = abundance ~ <15 unique rhs>, data = p.res.hab, 
##      family = negative.binomial(theta = 0.662012821216055), na.action = 
##      na.fail)
## 
## Component models: 
##        df   logLik    AICc delta weight
## 15      5 -1053.13 2116.58  0.00   0.16
## 135     6 -1052.12 2116.70  0.11   0.15
## 1235    7 -1051.29 2117.21  0.62   0.11
## 1345    7 -1051.62 2117.86  1.28   0.08
## 156     6 -1052.74 2117.95  1.36   0.08
## 145     6 -1053.00 2118.46  1.87   0.06
## 125     6 -1053.03 2118.53  1.94   0.06
## 1356    7 -1051.96 2118.54  1.96   0.06
## 12345   8 -1051.03 2118.86  2.28   0.05
## 12356   8 -1051.12 2119.05  2.46   0.05
## 13456   8 -1051.28 2119.38  2.79   0.04
## 1456    7 -1052.40 2119.43  2.84   0.04
## 1256    7 -1052.57 2119.77  3.19   0.03
## 123456  9 -1050.72 2120.46  3.88   0.02
## 1245    7 -1052.96 2120.54  3.95   0.02
## 
## Term codes: 
##  Habitat     s.do s.J.date    s.sal   s.temp     Year 
##        1        2        3        4        5        6 
## 
## Model-averaged coefficients:  
## (full average) 
##                  Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)       5.89148    0.21936     0.22079  26.684   <2e-16 ***
## HabitatSand flat -2.15828    0.34126     0.34339   6.285   <2e-16 ***
## s.temp            0.83582    0.13190     0.13263   6.302   <2e-16 ***
## s.J.date          0.10158    0.12777     0.12817   0.793    0.428    
## s.do              0.03780    0.08770     0.08803   0.429    0.668    
## s.sal            -0.04085    0.11108     0.11158   0.366    0.714    
## Year             -0.05453    0.14802     0.14871   0.367    0.714    
##  
## (conditional average) 
##                  Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)        5.8915     0.2194      0.2208  26.684   <2e-16 ***
## HabitatSand flat  -2.1583     0.3413      0.3434   6.285   <2e-16 ***
## s.temp             0.8358     0.1319      0.1326   6.302   <2e-16 ***
## s.J.date           0.1826     0.1206      0.1214   1.504    0.132    
## s.do               0.1103     0.1202      0.1209   0.912    0.362    
## s.sal             -0.1308     0.1666      0.1676   0.781    0.435    
## Year              -0.1749     0.2219      0.2233   0.783    0.434    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Relative variable importance: 
##                      Habitat s.temp s.J.date s.do s.sal Year
## Importance:          1.00    1.00   0.56     0.34 0.31  0.31
## N containing models:   15      15      8        7    7     7
p.res.ci <- data.frame(confint(p.avg_model.res.4, full = TRUE))

# Get pseudo R squared values for models up to delta < 4 model.list.res<-
# list(#manually list top x models from dredge)
p.res.4.Rsq <- rsquared(p.model.set.res.4, aicc = TRUE)

## write tables to .csv for easy comparison and plugging into appendix table
p.avg_model_4df.res <- data.frame(p.avg_model.res.4$msTable)
p.avg_model_components4.res <- cbind(p.res.4.Rsq, p.avg_model_4df.res)
p.r = data.frame(Coeff = rownames(p.avg_model_4df.res, rep(NA, length(p.avg_model_components4.res))))
p.avg_model_components4.res <- cbind(p.avg_model_components4.res, p.r)
p.avg_model_components4.res <- p.avg_model_components4.res[, -c(7, 8)]
# write.csv(p.avg_model_components4.res,
# '/Users/Lia/Documents/Git/Fraser-salmon/all.data/avg_model_components4_purseresidents.csv')

Parameter Plot

The results of model averaging including all top ranked models up to delta AICc 4

## List of 11
##  $ axis.title.x    : list()
##   ..- attr(*, "class")= chr [1:2] "element_blank" "element"
##  $ axis.title.y    :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : num 14
##   ..$ hjust        : NULL
##   ..$ vjust        : NULL
##   ..$ angle        : num 90
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi FALSE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text.x     : list()
##   ..- attr(*, "class")= chr [1:2] "element_blank" "element"
##  $ axis.text.y     :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : num 12
##   ..$ hjust        : NULL
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi FALSE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.ticks.x    : list()
##   ..- attr(*, "class")= chr [1:2] "element_blank" "element"
##  $ axis.line.x     :List of 6
##   ..$ colour       : chr "black"
##   ..$ size         : num 0.5
##   ..$ linetype     : NULL
##   ..$ lineend      : NULL
##   ..$ arrow        : logi FALSE
##   ..$ inherit.blank: logi FALSE
##   ..- attr(*, "class")= chr [1:2] "element_line" "element"
##  $ axis.line.y     :List of 6
##   ..$ colour       : chr "black"
##   ..$ size         : num 0.5
##   ..$ linetype     : NULL
##   ..$ lineend      : NULL
##   ..$ arrow        : logi FALSE
##   ..$ inherit.blank: logi FALSE
##   ..- attr(*, "class")= chr [1:2] "element_line" "element"
##  $ panel.border    : list()
##   ..- attr(*, "class")= chr [1:2] "element_blank" "element"
##  $ panel.grid.major: list()
##   ..- attr(*, "class")= chr [1:2] "element_blank" "element"
##  $ panel.grid.minor: list()
##   ..- attr(*, "class")= chr [1:2] "element_blank" "element"
##  $ plot.margin     :Classes 'margin', 'unit'  atomic [1:4] 2 2 2 2
##   .. ..- attr(*, "valid.unit")= int 8
##   .. ..- attr(*, "unit")= chr "pt"
##  - attr(*, "class")= chr [1:2] "theme" "gg"
##  - attr(*, "complete")= logi FALSE
##  - attr(*, "validate")= logi TRUE